Actual source code: ex29.c
1: /*$Id: $*/
3: /* solve the equations for the perturbed magnetic field only */
4: #define EQ
6: /* turning this on causes instability?!? */
7: #undef UPWINDING
9: static char help[] = "XXXXX with multigrid and timestepping 2d.n
10: n
11: -da_grid_x 6 -dmmg_nlevels 3 -da_grid_y 6 -mg_coarse_pc_type lu -mg_coarse_pc_lu_damping -mg_levels_pc_ilu_damping -mat_aij_no_inode n
12: -viscosity <nu>n
13: -skin_depth <d_e>n
14: -larmor_radius <rho_s>n
15: -contours : draw contour plots of solutionnn";
17: /*T
18: Concepts: SNES^solving a system of nonlinear equations (parallel multicomponent example);
19: Concepts: DA^using distributed arrays;
20: Concepts: multicomponent
21: Processors: n
22: T*/
24: /* ------------------------------------------------------------------------
26: We thank XXXXXX for contributing the FormFunctionLocal()
29: ------------------------------------------------------------------------- */
31: /*
32: Include "petscda.h" so that we can use distributed arrays (DAs).
33: Include "petscsnes.h" so that we can use SNES solvers.
34: Include "petscmg.h" to control the multigrid solvers.
35: Note that these automatically include:
36: petsc.h - base PETSc routines petscvec.h - vectors
37: petscsys.h - system routines petscmat.h - matrices
38: petscis.h - index sets petscksp.h - Krylov subspace methods
39: petscviewer.h - viewers petscpc.h - preconditioners
40: petscsles.h - linear solvers
41: */
42: #include petscsnes.h
43: #include petscda.h
44: #include petscmg.h
46: #define D_x(x,m,i,j) (p5 * (x[(j)][(i)+1].m - x[(j)][(i)+1].m) * dhx)
47: #define D_xm(x,m,i,j) ((x[(j)][(i)].m - x[(j)][(i)-1].m) * dhx)
48: #define D_xp(x,m,i,j) ((x[(j)][(i+1)].m - x[(j)][(i)].m) * dhx)
49: #define D_y(x,m,i,j) (p5 * (x[(j)+1][(i)].m - x[(j)-1][(i)].m) * dhy)
50: #define D_ym(x,m,i,j) ((x[(j)][(i)].m - x[(j)-1][(i)].m) * dhy)
51: #define D_yp(x,m,i,j) ((x[(j)+1][(i)].m - x[(j)][(i)].m) * dhy)
52: #define D_xx(x,m,i,j) ((x[(j)][(i)+1].m - two*x[(j)][(i)].m + x[(j)][(i)-1].m) * hydhx * dhxdhy)
53: #define D_yy(x,m,i,j) ((x[(j)+1][(i)].m - two*x[(j)][(i)].m + x[(j)-1][(i)].m) * hxdhy * dhxdhy)
54: #define Lapl(x,m,i,j) (D_xx(x,m,i,j) + D_yy(x,m,i,j))
55: #define lx (2.*M_PI)
56: #define ly (4.*M_PI)
57: #define sqr(a) ((a)*(a))
59: /*
60: User-defined routines and data structures
61: */
63: typedef struct {
64: PassiveScalar fnorm_ini,dt_ini;
65: PassiveScalar fnorm,dt,dt_out;
66: PassiveScalar ptime;
67: PassiveScalar max_time;
68: PassiveScalar fnorm_ratio;
69: int ires,itstep;
70: int max_steps,print_freq;
71: PassiveScalar t;
72: } TstepCtx;
74: typedef struct {
75: PetscScalar phi,U,psi,F;
76: } Field;
78: typedef struct {
79: PassiveScalar phi,U,psi,F;
80: } PassiveField;
82: typedef struct {
83: int mglevels;
84: int cycles; /* numbers of time steps for integration */
85: PassiveReal nu,d_e,rho_s; /* physical parameters */
86: PetscTruth draw_contours; /* flag - 1 indicates drawing contours */
87: PetscTruth PreLoading;
88: } Parameter;
90: typedef struct {
91: Vec Xold,func;
92: TstepCtx *tsCtx;
93: Parameter *param;
94: } AppCtx;
96: extern int FormFunctionLocal(DALocalInfo*,Field**,Field**,void*);
97: extern int Update(DMMG *);
98: extern int Initialize(DMMG *);
99: extern int AddTSTermLocal(DALocalInfo* info,Field **x,Field **f,AppCtx *user);
100: extern int Gnuplot(DA da, Vec X, double time);
101: extern int AttachNullSpace(PC,Vec);
103: #undef __FUNCT__
105: int main(int argc,char **argv)
106: {
107: DMMG *dmmg; /* multilevel grid structure */
108: AppCtx *user; /* user-defined work context (one for each level) */
109: TstepCtx tsCtx; /* time-step parameters (one total) */
110: Parameter param; /* physical parameters (one total) */
111: int i,ierr;
112: MPI_Comm comm;
113: DA da;
115: PetscInitialize(&argc,&argv,(char *)0,help);
116: comm = PETSC_COMM_WORLD;
119: PreLoadBegin(PETSC_TRUE,"SetUp");
121: param.PreLoading = PreLoading;
122: DMMGCreate(comm,1,&user,&dmmg);
123: param.mglevels = DMMGGetLevels(dmmg);
125: /*
126: Create distributed array multigrid object (DMMG) to manage parallel grid and vectors
127: for principal unknowns (x) and governing residuals (f)
128: */
129: DACreate2d(comm,DA_XYPERIODIC,DA_STENCIL_STAR,-6,-6,PETSC_DECIDE,PETSC_DECIDE,4,1,0,0,&da);
130: DMMGSetDM(dmmg,(DM)da);
131: DADestroy(da);
133: /*
134: Problem parameters
135: */
136: param.nu = 0.0;
137: param.rho_s = 0.0;
138: param.d_e = 0.2;
139: PetscOptionsGetReal(PETSC_NULL,"-viscosity",¶m.nu,PETSC_NULL);
140: PetscOptionsGetReal(PETSC_NULL,"-skin_depth",¶m.d_e,PETSC_NULL);
141: PetscOptionsGetReal(PETSC_NULL,"-larmor_radius",¶m.rho_s,PETSC_NULL);
142: PetscOptionsHasName(PETSC_NULL,"-contours",¶m.draw_contours);
144: DASetFieldName(DMMGGetDA(dmmg),0,"phi");
145: DASetFieldName(DMMGGetDA(dmmg),1,"U");
146: DASetFieldName(DMMGGetDA(dmmg),2,"psi");
147: DASetFieldName(DMMGGetDA(dmmg),3,"F");
149: /*======================================================================*/
150: /* Initilize stuff related to time stepping */
151: /*======================================================================*/
152: tsCtx.fnorm_ini = 0.0;
153: tsCtx.max_steps = 50;
154: tsCtx.max_time = 1.0e+12;
155: tsCtx.dt = .1;
156: tsCtx.fnorm_ratio = 1.0e+10;
157: tsCtx.t = 0.;
158: tsCtx.dt_out = .1;
159: PetscOptionsGetInt(PETSC_NULL,"-max_st",&tsCtx.max_steps,PETSC_NULL);
160: PetscOptionsGetReal(PETSC_NULL,"-ts_rtol",&tsCtx.fnorm_ratio,PETSC_NULL);
161: tsCtx.print_freq = tsCtx.max_steps;
162: PetscOptionsGetInt(PETSC_NULL,"-print_freq",&tsCtx.print_freq,PETSC_NULL);
163: PetscOptionsGetScalar(PETSC_NULL,"-deltat",&tsCtx.dt,PETSC_NULL);
165: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
166: Create user context, set problem data, create vector data structures.
167: Also, compute the initial guess.
168: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
169: /* create application context for each level */
170: PetscMalloc(param.mglevels*sizeof(AppCtx),&user);
171: for (i=0; i<param.mglevels; i++) {
172: /* create work vectors to hold previous time-step solution and function value */
173: VecDuplicate(dmmg[i]->x, &user[i].Xold);
174: VecDuplicate(dmmg[i]->x, &user[i].func);
175: user[i].tsCtx = &tsCtx;
176: user[i].param = ¶m;
177: dmmg[i]->user = &user[i];
178: }
179:
180: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
181: Create nonlinear solver context
182:
183: Process adiC: AddTSTermLocal FormFunctionLocal
184: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
185: DMMGSetSNESLocal(dmmg,FormFunctionLocal,0,ad_FormFunctionLocal,admf_FormFunctionLocal);
187: /* attach nullspace to each level of the preconditioner */
188: {
189: SLES subsles,sles;
190: PC pc,subpc;
191: PetscTruth mg;
193: SNESGetSLES(DMMGGetSNES(dmmg),&sles);
194: SLESGetPC(sles,&pc);
195: AttachNullSpace(pc,DMMGGetx(dmmg));
196: PetscTypeCompare((PetscObject)pc,PCMG,&mg);
197: if (mg) {
198: for (i=0; i<param.mglevels; i++) {
199: MGGetSmoother(pc,i,&subsles);
200: SLESGetPC(subsles,&subpc);
201: AttachNullSpace(subpc,dmmg[i]->x);
202: }
203: }
204: }
206: PetscPrintf(comm,"# viscosity = %g, skin_depth # = %g, larmor_radius # = %gn",
207: param.nu,param.d_e,param.rho_s);
208:
209:
210: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
211: Solve the nonlinear system
212: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
213:
214: PreLoadStage("Solve");
216: if (param.draw_contours) {
217: VecView(((AppCtx*)dmmg[param.mglevels-1]->user)->Xold,PETSC_VIEWER_DRAW_WORLD);
218: }
219: Update(dmmg);
221: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
222: Free work space. All PETSc objects should be destroyed when they
223: are no longer needed.
224: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
225:
226: for (i=0; i<param.mglevels; i++) {
227: VecDestroy(user[i].Xold);
228: VecDestroy(user[i].func);
229: }
230: PetscFree(user);
231: DMMGDestroy(dmmg);
232: PreLoadEnd();
233:
234: PetscFinalize();
235: return 0;
236: }
238: /* ------------------------------------------------------------------- */
239: #undef __FUNCT__
241: /* ------------------------------------------------------------------- */
242: int Gnuplot(DA da, Vec X, double time)
243: {
244: int i,j,xs,ys,xm,ym;
245: int xints,xinte,yints,yinte;
246: int ierr;
247: Field **x;
248: FILE *f;
249: char fname[100];
250: int cpu;
252: MPI_Comm_rank(PETSC_COMM_WORLD,&cpu);
253: sprintf(fname, "out-%g-%d.dat", time, cpu);
254: f = fopen(fname, "w");
255: DAGetCorners(da,&xs,&ys,PETSC_NULL,&xm,&ym,PETSC_NULL);
257: DAVecGetArray(da,X,(void**)&x);
259: xints = xs; xinte = xs+xm; yints = ys; yinte = ys+ym;
261: for (j=yints; j<yinte; j++) {
262: for (i=xints; i<xinte; i++) {
263: fprintf(f, "%d %d %g %g %g %g %g %gn", i, j, 0., 0., x[j][i].U, x[j][i].F, x[j][i].phi, x[j][i].psi);
264: }
265: fprintf(f, "n");
266: }
267: DAVecRestoreArray(da,X,(void**)&x);
268: fclose(f);
269: return 0;
270: }
272: /* ------------------------------------------------------------------- */
273: #undef __FUNCT__
275: /* ------------------------------------------------------------------- */
276: int Initialize(DMMG *dmmg)
277: {
278: AppCtx *appCtx = (AppCtx*)dmmg[0]->user;
279: Parameter *param = appCtx->param;
280: DA da;
281: int i,j,mx,my,ierr,xs,ys,xm,ym;
282: PetscReal two = 2.0,one = 1.0;
283: PetscReal hx,hy,dhx,dhy,hxdhy,hydhx,hxhy,dhxdhy;
284: PetscReal d_e,rho_s,de2,xx,yy;
285: Field **x, **localx;
286: Vec localX;
289: d_e = param->d_e;
290: rho_s = param->rho_s;
291: de2 = sqr(param->d_e);
293: da = (DA)(dmmg[param->mglevels-1]->dm);
294: DAGetInfo(da,0,&mx,&my,0,0,0,0,0,0,0,0);
296: dhx = mx/lx; dhy = my/ly;
297: hx = one/dhx; hy = one/dhy;
298: hxdhy = hx*dhy; hydhx = hy*dhx;
299: hxhy = hx*hy; dhxdhy = dhx*dhy;
301: /*
302: Get local grid boundaries (for 2-dimensional DA):
303: xs, ys - starting grid indices (no ghost points)
304: xm, ym - widths of local grid (no ghost points)
305: */
306: DAGetCorners(da,&xs,&ys,PETSC_NULL,&xm,&ym,PETSC_NULL);
308: DAGetLocalVector(da,&localX);
309: /*
310: Get a pointer to vector data.
311: - For default PETSc vectors, VecGetArray() returns a pointer to
312: the data array. Otherwise, the routine is implementation dependent.
313: - You MUST call VecRestoreArray() when you no longer need access to
314: the array.
315: */
316: DAVecGetArray(da,dmmg[param->mglevels-1]->x,(void**)&x);
317: DAVecGetArray(da,localX,(void**)&localx);
319: /*
320: Compute initial guess over the locally owned part of the grid
321: */
322: {
323: PetscReal eps = lx/ly;
324: PetscReal pert = 1e-4;
325: PetscReal k = 1.*eps;
326: PetscReal gam;
328: if (d_e < rho_s) d_e = rho_s;
329: gam = k * d_e;
331: for (j=ys-1; j<ys+ym+1; j++) {
332: yy = j * hy;
333: for (i=xs-1; i<xs+xm+1; i++) {
334: xx = i * hx;
336: if (xx < -M_PI/2) {
337: localx[j][i].phi = pert * gam / k * erf((xx + M_PI) / (sqrt(2) * d_e)) * (-sin(k*yy));
338: } else if (xx < M_PI/2) {
339: localx[j][i].phi = - pert * gam / k * erf(xx / (sqrt(2) * d_e)) * (-sin(k*yy));
340: } else if (xx < 3*M_PI/2){
341: localx[j][i].phi = pert * gam / k * erf((xx - M_PI) / (sqrt(2) * d_e)) * (-sin(k*yy));
342: } else {
343: localx[j][i].phi = - pert * gam / k * erf((xx - 2.*M_PI) / (sqrt(2) * d_e)) * (-sin(k*yy));
344: }
345: #ifdef EQ
346: localx[j][i].psi = 0.;
347: #else
348: localx[j][i].psi = cos(xx);
349: #endif
350: }
351: }
352: for (j=ys; j<ys+ym; j++) {
353: for (i=xs; i<xs+xm; i++) {
354: x[j][i].U = Lapl(localx,phi,i,j);
355: x[j][i].F = localx[j][i].psi - de2 * Lapl(localx,psi,i,j);
356: x[j][i].phi = localx[j][i].phi;
357: x[j][i].psi = localx[j][i].psi;
358: }
359: }
360: }
361: /*
362: Restore vector
363: */
364: DAVecRestoreArray(da,dmmg[param->mglevels-1]->x,(void**)&x);
365: DAVecRestoreArray(da,localX,(void**)&localx);
366: DARestoreLocalVector(da,&localX);
368: return(0);
369: }
371: #undef __FUNCT__
373: int ComputeMaxima(DA da, Vec X, PetscReal t)
374: {
375: int ierr,i,j,mx,my,xs,ys,xm,ym;
376: int xints,xinte,yints,yinte;
377: Field **x;
378: double norm[4];
379: double gnorm[4];
380: MPI_Comm comm;
383: DAGetInfo(da,0,&mx,&my,0,0,0,0,0,0,0,0);
385: DAGetCorners(da,&xs,&ys,PETSC_NULL,&xm,&ym,PETSC_NULL);
386: xints = xs; xinte = xs+xm; yints = ys; yinte = ys+ym;
388: DAVecGetArray(da,X,(void**)&x);
389: for (j=yints; j<yinte; j++) {
390: for (i=xints; i<xinte; i++) {
391: norm[0] = PetscMax(norm[0],x[j][i].U);
392: norm[1] = PetscMax(norm[1],x[j][i].F);
393: norm[2] = PetscMax(norm[2],x[j][i].phi);
394: norm[3] = PetscMax(norm[3],x[j][i].psi);
395: }
396: }
397: DAVecRestoreArray(da,X,(void**)&x);
398: PetscObjectGetComm((PetscObject)da,&comm);
399: MPI_Allreduce(norm, gnorm, 4, MPI_DOUBLE, MPI_MAX, comm);
400: fprintf(stderr, "%gt%gt%gt%gt%gn", t,
401: gnorm[0], gnorm[1], gnorm[2], gnorm[3]);
402: fflush(stderr);
404: return(0);
405: }
407: #undef __FUNCT__
409: int FormFunctionLocal(DALocalInfo *info,Field **x,Field **f,void *ptr)
410: {
411: AppCtx *user = (AppCtx*)ptr;
412: TstepCtx *tsCtx = user->tsCtx;
413: int ierr,i,j;
414: int xints,xinte,yints,yinte;
415: PassiveReal hx,hy,dhx,dhy,hxdhy,hydhx,hxhy,dhxdhy;
416: PassiveReal de2,rhos2,nu,dde2;
417: PassiveReal two = 2.0,one = 1.0,p5 = 0.5;
418: PetscScalar vx,vy,avx,avy,vxp,vxm,vyp,vym;
419: PetscScalar Bx,By,aBx,aBy,Bxp,Bxm,Byp,Bym;
420: PetscScalar xx;
421: PetscScalar F_eq_x;
422: PetscScalar By_eq;
425: de2 = sqr(user->param->d_e);
426: rhos2 = sqr(user->param->rho_s);
427: nu = user->param->nu;
428: dde2 = one/de2;
430: /*
431: Define mesh intervals ratios for uniform grid.
432: [Note: FD formulae below are normalized by multiplying through by
433: local volume element to obtain coefficients O(1) in two dimensions.]
434: */
435: dhx = info->mx/lx; dhy = info->my/ly;
436: hx = one/dhx; hy = one/dhy;
437: hxdhy = hx*dhy; hydhx = hy*dhx;
438: hxhy = hx*hy; dhxdhy = dhx*dhy;
440: xints = info->xs; xinte = info->xs+info->xm;
441: yints = info->ys; yinte = info->ys+info->ym;
443: /* Compute over the interior points */
444: for (j=yints; j<yinte; j++) {
445: for (i=xints; i<xinte; i++) {
446: #ifdef EQ
447: xx = i * hx;
448: F_eq_x = - (1. + de2) * sin(xx);
449: By_eq = sin(xx);
450: #else
451: F_eq_x = 0.;
452: By_eq = 0.;
453: #endif
455: /*
456: * convective coefficients for upwinding
457: */
459: vx = - D_y(x,phi,i,j);
460: vy = D_x(x,phi,i,j);
461: avx = PetscAbsScalar(vx); vxp = p5*(vx+avx); vxm = p5*(vx-avx);
462: avy = PetscAbsScalar(vy); vyp = p5*(vy+avy); vym = p5*(vy-avy);
463: #ifdef UPWINDING
464: vxp = vxm = p5*vx;
465: vyp = vym = p5*vy;
466: #endif
468: Bx = D_y(x,psi,i,j);
469: By = - D_x(x,psi,i,j) + By_eq;
470: aBx = PetscAbsScalar(Bx); Bxp = p5*(Bx+aBx); Bxm = p5*(Bx-aBx);
471: aBy = PetscAbsScalar(By); Byp = p5*(By+aBy); Bym = p5*(By-aBy);
472: #ifdef UPWINDING
473: Bxp = Bxm = p5*Bx;
474: Byp = Bym = p5*By;
475: #endif
477: /* Lap phi - U */
478: f[j][i].phi = (Lapl(x,phi,i,j) - x[j][i].U) * hxhy;
480: /* (1 - de^2 Lap) psi - F */
481: f[j][i].psi = (x[j][i].psi - de2 * Lapl(x,psi,i,j) - x[j][i].F) * hxhy;
483: /* - nu Lap U + vx * U_x + vy * U_y - (B_x F_x + B_y F_y)/de2 */
484: f[j][i].U = hxhy * (-nu * Lapl(x,U,i,j) +
485: (vxp*(D_xm(x,U,i,j)) +
486: vxm*(D_xp(x,U,i,j)) +
487: vyp*(D_ym(x,U,i,j)) +
488: vym*(D_yp(x,U,i,j))) -
489: (Bxp*(D_xm(x,F,i,j) + F_eq_x) +
490: Bxm*(D_xp(x,F,i,j) + F_eq_x) +
491: Byp*(D_ym(x,F,i,j)) +
492: Bym*(D_yp(x,F,i,j))) * dde2);
493:
494: /* -nu Lap F + vx * F_x + vy * F_y - rho_s2 (B_x U_x + B_y U_y) */
495: f[j][i].F = hxhy * (-nu * Lapl(x,F,i,j) + /* not quite right */
496: (vxp*(D_xm(x,F,i,j) + F_eq_x) +
497: vxm*(D_xp(x,F,i,j) + F_eq_x) +
498: vyp*(D_ym(x,F,i,j)) +
499: vym*(D_yp(x,F,i,j))) +
500: (Bxp*(D_xm(x,U,i,j)) +
501: Bxm*(D_xp(x,U,i,j)) +
502: Byp*(D_ym(x,U,i,j)) +
503: Bym*(D_yp(x,U,i,j))) * rhos2);
504: }
505: }
507: /* Add time step contribution */
508: if (tsCtx->ires) {
509: AddTSTermLocal(info,x,f,user);
510: }
511: /*
512: Flop count (multiply-adds are counted as 2 operations)
513: */
514: /* PetscLogFlops(84*info->ym*info->xm); FIXME */
515: return(0);
516: }
518: /*---------------------------------------------------------------------*/
519: #undef __FUNCT__
521: int Update(DMMG *dmmg)
522: /*---------------------------------------------------------------------*/
523: {
524:
525: AppCtx *user = (AppCtx *) ((dmmg[0])->user);
526: TstepCtx *tsCtx = user->tsCtx;
527: Parameter *param = user->param;
528: SNES snes;
529: int ierr,its,lits,i;
530: int max_steps;
531: PetscTruth print_flag = PETSC_FALSE;
532: int nfailsCum = 0,nfails = 0;
533: static int ic_out;
537: PetscOptionsHasName(PETSC_NULL,"-print",&print_flag);
538: if (user->param->PreLoading)
539: max_steps = 1;
540: else
541: max_steps = tsCtx->max_steps;
542:
543: Initialize(dmmg);
545: for (tsCtx->itstep = 0; tsCtx->itstep < max_steps; tsCtx->itstep++) {
546: for (i=param->mglevels-1; i>0 ;i--) {
547: MatRestrict(dmmg[i]->R, dmmg[i]->x, dmmg[i-1]->x);
548: VecPointwiseMult(dmmg[i]->Rscale,dmmg[i-1]->x,dmmg[i-1]->x);
549: VecCopy(dmmg[i]->x, ((AppCtx*)dmmg[i]->user)->Xold);
550: }
551: VecCopy(dmmg[0]->x, ((AppCtx*)dmmg[0]->user)->Xold);
553: DMMGSolve(dmmg);
554: snes = DMMGGetSNES(dmmg);
555: SNESGetIterationNumber(snes,&its);
556: SNESGetNumberLinearIterations(snes,&lits);
557: PetscPrintf(PETSC_COMM_WORLD,"Number of Newton iterations = %d / lin = %dn", its,lits);
558: SNESGetNumberUnsuccessfulSteps(snes,&nfails);
559: PetscPrintf(PETSC_COMM_WORLD,"Number of unsuccessful = %dn", nfails);
560: nfailsCum += nfails; nfails = 0;
561: if (nfailsCum >= 2) SETERRQ(1,"Unable to find a Newton Step");
563: SNESComputeFunction(snes,dmmg[param->mglevels-1]->x,((AppCtx*)dmmg[param->mglevels-1]->user)->func);
564: VecNorm(user->func,NORM_2,&tsCtx->fnorm);
565:
566: tsCtx->t += tsCtx->dt;
567: if (print_flag) {
568: PetscPrintf(PETSC_COMM_WORLD,"After Time Step %d and fnorm = %gn",
569: tsCtx->itstep,tsCtx->fnorm);
570: }
571: if (!param->PreLoading) {
572: if (param->draw_contours) {
573: VecView(DMMGGetx(dmmg),PETSC_VIEWER_DRAW_WORLD);
574: }
575: /* compute maxima */
576: ComputeMaxima((DA) dmmg[param->mglevels-1]->dm, dmmg[param->mglevels-1]->x, tsCtx->t);
577: /* output */
578: if (ic_out++ == (int)(tsCtx->dt_out / tsCtx->dt + .5)) {
579: Gnuplot((DA) dmmg[param->mglevels-1]->dm,
580: dmmg[param->mglevels-1]->x, tsCtx->t);
581: ic_out = 1;
582: }
583: }
584: } /* End of time step loop */
585:
586: PetscPrintf(PETSC_COMM_WORLD,"timesteps %d fnorm = %gn",tsCtx->itstep,tsCtx->fnorm);
587: if (user->param->PreLoading) {
588: tsCtx->fnorm_ini = 0.0;
589: }
590:
591: return(0);
592: }
594: /*---------------------------------------------------------------------*/
595: #undef __FUNCT__
597: int AddTSTermLocal(DALocalInfo* info,Field **x,Field **f,AppCtx *user)
598: /*---------------------------------------------------------------------*/
599: {
600: TstepCtx *tsCtx = user->tsCtx;
601: DA da = info->da;
602: int ierr,i,j;
603: int xints,xinte,yints,yinte;
604: PassiveReal hx,hy,dhx,dhy,hxhy;
605: PassiveReal one = 1.0;
606: PassiveScalar dtinv;
607: PassiveField **xold;
611: xints = info->xs; xinte = info->xs+info->xm;
612: yints = info->ys; yinte = info->ys+info->ym;
614: dhx = info->mx/lx; dhy = info->my/ly;
615: hx = one/dhx; hy = one/dhy;
616: hxhy = hx*hy;
618: ierr = DAVecGetArray(da,user->Xold,(void**)&xold);
619: dtinv = hxhy/(tsCtx->dt);
620: for (j=yints; j<yinte; j++) {
621: for (i=xints; i<xinte; i++) {
622: f[j][i].U += dtinv*(x[j][i].U-xold[j][i].U);
623: f[j][i].F += dtinv*(x[j][i].F-xold[j][i].F);
624: }
625: }
626: DAVecRestoreArray(da,user->Xold,(void**)&xold);
627: return(0);
628: }
630: /* -------------------------------------------------------------------------*/
631: #undef __FUNCT__
633: int AttachNullSpace(PC pc,Vec model)
634: {
635: int i,ierr,rstart,rend,N;
636: MatNullSpace sp;
637: Vec v,vs[1];
638: PetscScalar *vx,scale;
641: ierr = VecDuplicate(model,&v);
642: ierr = VecGetSize(model,&N);
643: scale = 2.0/sqrt(N);
644: ierr = VecGetArray(v,&vx);
645: ierr = VecGetOwnershipRange(v,&rstart,&rend);
646: for (i=rstart; i<rend; i++) {
647: if (!(i % 4)) vx[i-rstart] = scale;
648: else vx[i-rstart] = 0.0;
649: }
650: ierr = VecRestoreArray(v,&vx);
651: vs[0] = v;
652: ierr = MatNullSpaceCreate(PETSC_COMM_WORLD,0,1,vs,&sp);
653: ierr = VecDestroy(v);
654: ierr = PCNullSpaceAttach(pc,sp);
655: ierr = MatNullSpaceDestroy(sp);
656: return(0);
657: }