Actual source code: ex1.c
1: static char help[] = "\n";
2: /* - - - - - - - - - - - - - - - - - - - - - - - -
3: ex1.c
4: Simple example with 3 dof da: u,w,phi
5: where u & w are time-independent & analytically prescribed.
6: phi is advected explicitly.
7: - - - - - - - - - - - - - - - - - - - - - - - - */
8: #include <petscsnes.h>
9: #include <petscdmda.h>
10: #include <petscdmmg.h>
11: #include <petscbag.h>
12: #include <petsccharacteristic.h>
14: #define SHEAR_CELL 0
15: #define SOLID_BODY 1
16: #define FNAME_LENGTH 60
18: typedef struct field_s {
19: PetscReal u,w,phi;
20: } Field;
22: typedef struct parameter_s {
23: int ni, nj, pi, pj;
24: PetscReal sigma,xctr,zctr,L1,L2,LINF;
25: int verify_result, flow_type, sl_event;
26: PetscBool verify, param_test, output_to_file;
27: char output_filename[FNAME_LENGTH];
28: /* timestep stuff */
29: PetscReal t; /* the time */
30: int n; /* the time step */
31: PetscReal t_max, dt, cfl, t_output_interval;
32: int N_steps;
33: } Parameter;
35: typedef struct gridinfo_s {
36: DMDABoundaryType periodic;
37: DMDAStencilType stencil;
38: int ni,nj,dof,stencil_width,mglevels;
39: PetscReal dx,dz;
40: } GridInfo;
42: typedef struct appctx_s {
43: Vec Xold;
44: PetscBag bag;
45: GridInfo *grid;
46: } AppCtx;
48: /* Main routines */
49: int SetParams (AppCtx*);
50: int ReportParams (AppCtx*);
51: int Initialize (DMMG*);
52: int DoSolve (DMMG*);
53: int DoOutput (DMMG*, int);
54: int CalcSolnNorms (DMMG*, PetscReal*);
55: int DoVerification (DMMG*, AppCtx*);
56: int DMDASetFieldNames (const char*, const char*, const char*, DM);
57: PetscReal BiCubicInterp (Field**, PetscReal, PetscReal);
58: PetscReal CubicInterp (PetscReal, PetscReal, PetscReal, PetscReal, PetscReal);
59: PetscBool OptionsHasName(const char*);
61: /* characteristic call-backs (static interface) */
62: PetscErrorCode InterpVelocity2D(void*, PetscReal[], PetscInt, PetscInt[], PetscReal[], void*);
63: PetscErrorCode InterpFields2D (void*, PetscReal[], PetscInt, PetscInt[], PetscReal[], void*);
65: /* a few macros for convenience */
66: #define REG_REAL(A,B,C,D,E) ierr=PetscBagRegisterReal(A,B,C,D,E);CHKERRQ(ierr)
67: #define REG_INTG(A,B,C,D,E) ierr=PetscBagRegisterInt(A,B,C,D,E);CHKERRQ(ierr)
68: #define REG_TRUE(A,B,C,D,E) ierr=PetscBagRegisterBool(A,B,C,D,E);CHKERRQ(ierr)
69: #define REG_STRG(A,B,C,D,E,F) ierr=PetscBagRegisterString(A,B,C,D,E,F);CHKERRQ(ierr)
70: #define REG_ENUM(A,B,C,D,E,F) ierr=PetscBagRegisterEnum(A,B,C,D,E,F);CHKERRQ(ierr)
72: /*-----------------------------------------------------------------------*/
75: int main(int argc,char **argv)
76: /*-----------------------------------------------------------------------*/
77: {
78: DMMG *dmmg; /* multilevel grid structure */
79: AppCtx *user; /* user-defined work context */
80: Parameter *param;
81: GridInfo grid;
82: int ierr,result = 0;
83: MPI_Comm comm;
84: DM da;
86: PetscInitialize(&argc,&argv,(char *)0,help);
87: comm = PETSC_COMM_WORLD;
89: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
90: Set up the problem parameters.
91: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
92: PetscMalloc(sizeof(AppCtx),&user);
93: PetscBagCreate(comm,sizeof(Parameter),&(user->bag));
94: user->grid = &grid;
95: SetParams(user);
96: ReportParams(user);
97: PetscBagGetData(user->bag,(void**)¶m);
99: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
100: Create distributed array multigrid object (DMMG) to manage parallel grid and vectors
101: for principal unknowns (x) and governing residuals (f)
102: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
103: DMMGCreate(comm,grid.mglevels,user,&dmmg);
104: DMDACreate2d(comm,grid.periodic,grid.periodic,grid.stencil,grid.ni,grid.nj,PETSC_DECIDE,PETSC_DECIDE,grid.dof,grid.stencil_width,0,0,&da);
105: DMMGSetDM(dmmg,(DM)da);
106: DMDAGetInfo(da,PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL,&(param->pi),&(param->pj),PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL);
107: REG_INTG(user->bag,¶m->pi,param->pi ,"procs_x","<DO NOT SET> Processors in the x-direction");
108: REG_INTG(user->bag,¶m->pj,param->pj ,"procs_y","<DO NOT SET> Processors in the y-direction");
110: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
111: Create user context, set problem data, create vector data structures.
112: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
113: DMGetGlobalVector(da, &(user->Xold));
115: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
116: Initialize and solve the nonlinear system
117: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
118: Initialize(dmmg);
119: DoSolve(dmmg);
120: if (param->verify) result = param->verify_result;
121:
122: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
123: Free work space.
124: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
125: DMRestoreGlobalVector(da, &(user->Xold));
126: PetscBagDestroy(&user->bag);
127: PetscFree(user);
128: DMDestroy(&da);
129: DMMGDestroy(dmmg);
130: PetscFinalize();
131: return result;
132: }
134: /*---------------------------------------------------------------------*/
137: int SetParams(AppCtx *user)
138: /*---------------------------------------------------------------------*/
139: {
140: PetscBag bag = user->bag;
141: Parameter *p;
142: GridInfo *grid = user->grid;
143: int ierr, ierr_out=0;
144: PetscBagGetData(bag,(void**)&p);
146: /* give the bag a name */
147: PetscBagSetName(bag,"ex1_params","Parameter bag for ex1.c");
149: /* verification */
150: REG_TRUE(bag,&p->verify,PETSC_FALSE ,"verify","Do verification run (T/F)");
152: /* domain geometry & grid size */
153: REG_INTG(bag,&p->ni,40 ,"ni","Grid points in x-dir");
154: REG_INTG(bag,&p->nj,40 ,"nj","Grid points in y-dir");
155: grid->dx = 1.0/((double)(p->ni - 1));
156: grid->dz = 1.0/((double)(p->nj - 1));
158: /* initial conditions */
159: REG_INTG(bag,&p->flow_type,SHEAR_CELL ,"flow_type","Flow field mode: 0=shear cell, 1=translation");
160: REG_REAL(bag,&p->sigma,0.07 ,"sigma","Standard deviation of the gaussian IC");
161: REG_REAL(bag,&p->xctr,0.5 ,"xctr","x-position of the center of the gaussian IC");
162: REG_REAL(bag,&p->zctr,0.75 ,"zctr","z-position of the center of the gaussian IC");
164: /* time stepping */
165: REG_REAL(bag,&p->t_max,1 ,"t_max","Maximum dimensionless time");
166: REG_REAL(bag,&p->cfl,5 ,"cfl","Courant number");
167: REG_REAL(bag,&p->t_output_interval,0.1 ,"t_output","Dimensionless time interval for output");
168: REG_INTG(bag,&p->N_steps,1000 ,"nsteps","Maximum time-steps");
169: REG_INTG(bag,&p->n,1 ,"nsteps","<DO NOT SET> current time-step");
170: REG_REAL(bag,&p->t,0.0 ,"time","<DO NOT SET> initial time");
171: REG_REAL(bag,&p->dt,p->cfl*PetscMin(grid->dx,grid->dz),"dt","<DO NOT SET> time-step size");
173: /* output options */
174: REG_TRUE(bag,&p->param_test ,PETSC_FALSE ,"test","Run parameter test only (T/F)");
175: REG_STRG(bag,&p->output_filename,FNAME_LENGTH ,"null","output_file","Name base for output files, set with: -output_file <filename>");
176: REG_TRUE(bag,&p->output_to_file,PETSC_FALSE ,"do_output","<DO NOT SET> flag will be true if you specify an output file name");
177: p->output_to_file = OptionsHasName("-output_file");
179: if (p->verify) {
180: REG_INTG(bag,&p->verify_result,0 ,"ver_result","<DO NOT SET> Result of verification test");
181: REG_REAL(bag,&p->L1 ,4924.42 ,"L1","<DO NOT SET> L1");
182: REG_REAL(bag,&p->L2 ,496.287 ,"L2","<DO NOT SET> L2");
183: REG_REAL(bag,&p->LINF,100 ,"L3","<DO NOT SET> L3");
184:
185: p->verify_result = 0; p->L1 = 4924.42; p->L2 = 496.287; p->LINF = 100;
186: grid->ni = grid->nj = 40; grid->dx = 1.0/((double)(grid->ni)); grid->dz = 1.0/((double)(grid->nj));
187: p->flow_type = SHEAR_CELL; p->sigma = 0.07; p->xctr = 0.5; p->zctr = 0.75;
188: p->t_max = 0.5; p->cfl = 5; p->t_output_interval = 0.1; p->dt = p->cfl*PetscMin(grid->dx,grid->dz);
189: }
191: grid->ni = p->ni;
192: grid->nj = p->nj;
193: grid->periodic = DMDA_BOUNDARY_PERIODIC;
194: grid->stencil = DMDA_STENCIL_BOX;
195: grid->dof = 3;
196: grid->stencil_width = 2;
197: grid->mglevels = 1;
198: return ierr_out;
199: }
201: /*---------------------------------------------------------------------*/
204: int ReportParams(AppCtx *user)
205: /*---------------------------------------------------------------------*/
206: {
207: Parameter *param;
208: GridInfo *grid = user->grid;
209: int ierr, ierr_out=0;
210: PetscBagGetData(user->bag,(void**)¶m);
212: PetscPrintf(PETSC_COMM_WORLD,"---------------MOC test 1----------------\n");
213: PetscPrintf(PETSC_COMM_WORLD,"Prescribed wind, method of\n");
214: PetscPrintf(PETSC_COMM_WORLD,"characteristics advection, explicit time-\n");
215: PetscPrintf(PETSC_COMM_WORLD,"stepping.\n\n");
216: if (param->flow_type == 0) {
217: PetscPrintf(PETSC_COMM_WORLD,"Flow_type: %d (shear cell).\n\n",param->flow_type);
218: }
219: if (param->flow_type == 1) {
220: PetscPrintf(PETSC_COMM_WORLD,"Flow_type: %d (rigid body rotation).\n\n",param->flow_type);
221: }
222: if (param->verify) {
223: PetscPrintf(PETSC_COMM_WORLD," ** VERIFICATION RUN ** \n\n");
224: }
225: PetscPrintf(PETSC_COMM_WORLD," [ni,nj] = %d, %d [dx,dz] = %5.4g, %5.4g\n",grid->ni,grid->nj,grid->dx,grid->dz);
226: PetscPrintf(PETSC_COMM_WORLD," t_max = %g, cfl = %g, dt = %5.4g,",param->t_max,param->cfl,param->dt);
227: PetscPrintf(PETSC_COMM_WORLD," t_output = %g\n",param->t_output_interval);
228: if (param->output_to_file) {
229: PetscPrintf(PETSC_COMM_WORLD,"Output File: Binary file \"%s\"\n",param->output_filename);
230: }
231: if (!param->output_to_file)
232: PetscPrintf(PETSC_COMM_WORLD,"Output File: NO OUTPUT!\n");
233: PetscPrintf(PETSC_COMM_WORLD,"----------------------------------------\n");
234: if (param->param_test) PetscEnd();
235: return ierr_out;
236: }
238: /* ------------------------------------------------------------------- */
241: int Initialize(DMMG *dmmg)
242: /* ------------------------------------------------------------------- */
243: {
244: AppCtx *user = (AppCtx*)dmmg[0]->user;
245: Parameter *param;
246: DM da;
247: PetscReal PI = 3.14159265358979323846;
248: PetscReal sigma,xc,zc;
249: PetscReal dx=user->grid->dx,dz=user->grid->dz;
250: int i,j,ierr,is,js,im,jm;
251: Field **x;
252: PetscBagGetData(user->bag,(void**)¶m);
253: sigma=param->sigma; xc=param->xctr; zc=param->zctr;
255: /* Get the DM and grid */
256: da = (dmmg[0]->dm);
257: DMDAGetCorners(da,&is,&js,PETSC_NULL,&im,&jm,PETSC_NULL);
258: DMDAVecGetArray(da,user->Xold,(void**)&x);
260: for (j=js; j<js+jm; j++) {
261: for (i=is; i<is+im; i++) {
262: if (param->flow_type == SHEAR_CELL) {
263: x[j][i].u = -sin(PI*i*dx)*cos(PI*j*dz)/dx;
264: x[j][i].w = sin(PI*j*dz)*cos(PI*i*dx)/dz;
265: } else {
266: x[j][i].u = 0.0;
267: x[j][i].w = -1.0/dz;
268: }
269: x[j][i].phi = 100*exp(-0.5*((i*dx-xc)*(i*dx-xc)+(j*dz-zc)*(j*dz-zc))/sigma/sigma);
270: }
271: }
272:
273: /* restore the grid to it's vector */
274: DMDAVecRestoreArray(da,user->Xold,(void**)&x);
275: VecCopy(user->Xold, DMMGGetx(dmmg));
276: return 0;
277: }
279: /* ------------------------------------------------------------------- */
282: int DoSolve(DMMG *dmmg)
283: /* ------------------------------------------------------------------- */
284: {
285: AppCtx *user = (AppCtx*)dmmg[0]->user;
286: Parameter *param;
287: PetscReal t_output = 0.0;
288: int ierr, n_plot = 0, Ncomponents, components[3];
289: DM da = DMMGGetDM(dmmg);
290: Vec Xstar;
291: Characteristic c;
292: PetscBagGetData(user->bag,(void**)¶m);
294: DMGetGlobalVector(da, &Xstar);
296: /*------------ BEGIN CHARACTERISTIC SETUP ---------------*/
297: CharacteristicCreate(PETSC_COMM_WORLD, &c);
298: /* set up the velocity interpolation system */
299: Ncomponents = 2; components[0] = 0; components[1] = 1;
300: CharacteristicSetVelocityInterpolationLocal(c, da, DMMGGetx(dmmg), user->Xold, Ncomponents, components, InterpVelocity2D, user);
301: /* set up the fields interpolation system */
302: Ncomponents = 1; components[0] = 2;
303: CharacteristicSetFieldInterpolationLocal(c, da, user->Xold, Ncomponents, components, InterpFields2D, user);
304: /*------------ END CHARACTERISTIC SETUP ----------------*/
306: /* output initial data */
307: PetscPrintf(PETSC_COMM_WORLD," Initialization, Time: %5.4g\n", param->t);
308: if (param->verify) { DoVerification(dmmg,user); }
309: DoOutput(dmmg,n_plot);
310: t_output += param->t_output_interval; n_plot++;
312: /* timestep loop */
313: for (param->t=param->dt; param->t<=param->t_max; param->t+=param->dt) {
314: if (param->n > param->N_steps) {
315: PetscPrintf(PETSC_COMM_WORLD,"EXCEEDED MAX NUMBER OF TIMESTEPS! EXITING SOLVE!\n");
316: return 0;
317: }
319: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
320: Solve at time t & copy solution into solution vector.
321: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
322: /* Copy in the velocities to Xstar */
323: VecCopy(DMMGGetx(dmmg), Xstar);
324: /* Put \phi_* into Xstar */
325: CharacteristicSolve(c, param->dt, Xstar);
326: /* Copy the advected field into the solution \phi_t = \phi_* */
327: VecCopy(Xstar, DMMGGetx(dmmg));
329: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
330: Copy new solution to old solution in prep for the next timestep.
331: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
332: VecCopy(DMMGGetx(dmmg), user->Xold);
335: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
336: Timestep complete, report and update counter.
337: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
338: PetscPrintf(PETSC_COMM_WORLD," Step: %d, Time: %5.4g\n", param->n, param->t);
339: param->n++;
341: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
342: Verify and make output.
343: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
344: if (param->verify) { DoVerification(dmmg,user); }
345: if (param->t >= t_output) {
346: DoOutput(dmmg,n_plot);
347: t_output += param->t_output_interval; n_plot++;
348: }
349: }
350: DMRestoreGlobalVector(da, &Xstar);
351: CharacteristicDestroy(&c);
352: return 0;
353: }
355: /*---------------------------------------------------------------------*/
358: /* uses analytic velocity fields */
359: PetscErrorCode InterpVelocity2D(void *f, PetscReal ij_real[], PetscInt numComp,
360: PetscInt components[], PetscReal velocity[],
361: void *ctx)
362: /*---------------------------------------------------------------------*/
363: {
364: AppCtx *user = (AppCtx *) ctx;
365: Parameter *param;
366: PetscReal dx=user->grid->dx, dz=user->grid->dz;
367: PetscReal PI = 3.14159265358979323846;
368: int ierr;
369: PetscBagGetData(user->bag,(void**)¶m);
371: if (param->flow_type == SHEAR_CELL) {
372: velocity[0] = -sin(PI*ij_real[0]*dx)*cos(PI*ij_real[1]*dz)/dx;
373: velocity[1] = sin(PI*ij_real[1]*dz)*cos(PI*ij_real[0]*dx)/dz;
374: } else {
375: velocity[0] = 0.;
376: velocity[1] = -1./dz;
377: }
378: return 0;
379: }
381: /*---------------------------------------------------------------------*/
384: /* uses bicubic interpolation */
385: PetscErrorCode InterpFields2D(void *f, PetscReal ij_real[], PetscInt numComp,
386: PetscInt components[], PetscReal field[],
387: void *ctx)
388: /*---------------------------------------------------------------------*/
389: {
390: AppCtx *user = (AppCtx*)ctx;
391: Field **x = (Field**)f;
392: int ni=user->grid->ni, nj=user->grid->nj;
393: PetscReal ir=ij_real[0], jr=ij_real[1];
395: /* boundary condition: set to zero if characteristic begins outside the domain */
396: if ( ir < 0 || ir > ni-1 || jr < 0 || jr> nj-1 ) {
397: field[0] = 0.0;
398: } else {
399: field[0] = BiCubicInterp(x, ir, jr);
400: }
401: return 0;
402: }
404: /*---------------------------------------------------------------------*/
407: PetscReal BiCubicInterp(Field **x, PetscReal ir, PetscReal jr)
408: /*---------------------------------------------------------------------*/
409: {
410: int im, jm, imm,jmm,ip,jp,ipp,jpp;
411: PetscReal il, jl, row1, row2, row3, row4;
412: im = (int)floor(ir); jm = (int)floor(jr);
413: il = ir - im + 1.0; jl = jr - jm + 1.0;
414: imm = im-1; ip = im+1; ipp = im+2;
415: jmm = jm-1; jp = jm+1; jpp = jm+2;
416: row1 = CubicInterp(il,x[jmm][imm].phi,x[jmm][im].phi,x[jmm][ip].phi,x[jmm][ipp].phi);
417: row2 = CubicInterp(il,x[jm] [imm].phi,x[jm] [im].phi,x[jm] [ip].phi,x[jm] [ipp].phi);
418: row3 = CubicInterp(il,x[jp] [imm].phi,x[jp] [im].phi,x[jp] [ip].phi,x[jp] [ipp].phi);
419: row4 = CubicInterp(il,x[jpp][imm].phi,x[jpp][im].phi,x[jpp][ip].phi,x[jpp][ipp].phi);
420: return CubicInterp(jl,row1,row2,row3,row4);
421: }
423: /*---------------------------------------------------------------------*/
426: PetscReal CubicInterp(PetscReal x, PetscReal y_1, PetscReal y_2,
427: PetscReal y_3, PetscReal y_4)
428: /*---------------------------------------------------------------------*/
429: {
430: PetscReal sxth=0.16666666666667, retval;
431: retval = - y_1*(x-1.0)*(x-2.0)*(x-3.0)*sxth + y_2*(x-0.0)*(x-2.0)*(x-3.0)*0.5
432: - y_3*(x-0.0)*(x-1.0)*(x-3.0)*0.5 + y_4*(x-0.0)*(x-1.0)*(x-2.0)*sxth;
433: return retval;
434: }
436: /*---------------------------------------------------------------------*/
439: int DoVerification(DMMG *dmmg, AppCtx *user)
440: /*---------------------------------------------------------------------*/
441: {
442: Parameter *param;
443: PetscReal t1,t2,t3,norms[3];
444: int ierr;
445: PetscBagGetData(user->bag,(void**)¶m);
446:
447: CalcSolnNorms(dmmg, norms);
448: t1 = (norms[0]-param->L1)/param->L1*100.0;
449: t2 = (norms[1]-param->L2)/param->L2*100.0;
450: t3 = (norms[2]-param->LINF)/param->LINF*100.0;
451: if ((fabs(t1)>1.0) || (fabs(t2)>1.0) || (fabs(t3)>5.0)) {
452: param->verify_result = 1;
453: }
454: PetscPrintf(PETSC_COMM_WORLD," Step: %d, Soln norms: %g (L1) %g (L2) %g (LINF)\n",param->n,norms[0],norms[1],norms[2]);
455: PetscPrintf(PETSC_COMM_WORLD," Step: %d, Soln norms %%err: %5.2g (L1) %5.2g (L2) %5.2g (LINF)\n",param->n,t1,t2,t3);
456: return 0;
457: }
459: /*---------------------------------------------------------------------*/
462: int CalcSolnNorms(DMMG *dmmg, PetscReal norms[])
463: /*---------------------------------------------------------------------*/
464: {
465: Vec x;
466: int ierr;
467: x = DMMGGetx(dmmg);
468: VecNorm(x, NORM_1, &(norms[0]));
469: VecNorm(x, NORM_2, &(norms[1]));
470: VecNorm(x, NORM_INFINITY, &(norms[2]));
471: return 0;
472: }
474: /*---------------------------------------------------------------------*/
477: int DoOutput(DMMG *dmmg, int n_plot)
478: /*---------------------------------------------------------------------*/
479: {
480: AppCtx *user = (AppCtx*)dmmg[0]->user;
481: Parameter *param;
482: int ierr;
483: char filename[FNAME_LENGTH];
484: PetscViewer viewer;
485: DM da;
486: PetscBagGetData(user->bag,(void**)¶m);
487: da = DMMGGetDM(dmmg);
489: if (param->output_to_file) { /* send output to binary file */
490: /* generate filename for time t */
491: sprintf(filename,"%s_%3.3d",param->output_filename,n_plot);
492: PetscPrintf(PETSC_COMM_WORLD,"Generating output: time t = %g, ",param->t);
493: PetscPrintf(PETSC_COMM_WORLD,"file = \"%s\"\n",filename);
495: /* make output files */
496: PetscViewerBinaryMatlabOpen(PETSC_COMM_WORLD,filename,&viewer);
497: PetscViewerBinaryMatlabOutputBag(viewer,"par",user->bag);
498: DMDASetFieldNames("u","v","phi",da);
499: PetscViewerBinaryMatlabOutputVecDA(viewer,"field",DMMGGetx(dmmg),da);
500: PetscViewerBinaryMatlabDestroy(&viewer);
501: }
502: return 0;
503: }
505: /* ------------------------------------------------------------------- */
508: int DMDASetFieldNames(const char n0[], const char n1[], const char n2[], DM da)
509: /* ------------------------------------------------------------------- */
510: {
511: int ierr;
512: DMDASetFieldName(da,0,n0);
513: DMDASetFieldName(da,1,n1);
514: DMDASetFieldName(da,2,n2);
515: return 0;
516: }
518: /* ------------------------------------------------------------------- */
521: PetscBool OptionsHasName(const char name[])
522: /* ------------------------------------------------------------------- */
523: {
524: PetscBool retval;
525: int ierr;
526: PetscOptionsHasName(PETSC_NULL,name,&retval);
527: return retval;
528: }