Actual source code: ex2.c
1: static char help[] = "\n";
2: /* - - - - - - - - - - - - - - - - - - - - - - - -
3: ex2.c
4: Explicit semi-Lagrangian advection of a set of
5: passive tracers. DM has 7 dof: u,w,c1,c2,c3,c4
6: where u & w are time-independent & analytically prescribed.
7: This is identical to ex1 but with more advected
8: fields.
9: - - - - - - - - - - - - - - - - - - - - - - - - */
10: #include <petscsnes.h>
11: #include <petscdmda.h>
12: #include <petscdmmg.h>
13: #include <petscbag.h>
14: #include <petsccharacteristic.h>
16: #define SHEAR_CELL 0
17: #define SOLID_BODY 1
18: #define FNAME_LENGTH 60
20: typedef struct field_s {
21: PetscReal u,w,c1,c2,c3,c4;
22: } Field;
24: typedef PetscReal FieldArr[6];
26: typedef struct parameter_s {
27: int ni, nj, pi, pj;
28: PetscReal sigma,xctr,zctr,L1,L2,LINF;
29: int flow_type;
30: PetscBool param_test, output_to_file;
31: char output_filename[FNAME_LENGTH];
32: /* timestep stuff */
33: PetscReal t; /* the time */
34: int n; /* the time step */
35: PetscReal t_max, dt, cfl, t_output_interval;
36: int N_steps;
37: } Parameter;
39: typedef struct gridinfo_s {
40: DMDABoundaryType periodic;
41: DMDAStencilType stencil;
42: int ni,nj,dof,stencil_width,mglevels;
43: PetscReal dx,dz;
44: } GridInfo;
46: typedef struct appctx_s {
47: Vec Xold;
48: PetscBag bag;
49: GridInfo grid;
50: } AppCtx;
52: /* Main routines */
53: int SetParams (AppCtx*);
54: int ReportParams (AppCtx*);
55: int Initialize (DMMG*);
56: int DoSolve (DMMG*);
57: int DoOutput (DMMG*, int);
58: int DMDASetFieldNames (const char*,const char*,const char*,const char*,const char*,const char*,DM);
59: PetscReal BiCubicInterp (FieldArr**, PetscReal, PetscReal, int);
60: PetscReal CubicInterp (PetscReal, PetscReal, PetscReal, PetscReal, PetscReal);
61: PetscBool OptionsHasName(const char*);
63: /* characteristic call-backs (static interface) */
64: PetscErrorCode InterpVelocity2D(void*, PetscReal[], PetscInt, PetscInt[], PetscReal[], void*);
65: PetscErrorCode InterpFields2D (void*, PetscReal[], PetscInt, PetscInt[], PetscReal[], void*);
67: /* a few macros for convenience */
68: #define REG_REAL(A,B,C,D,E) ierr=PetscBagRegisterReal(A,B,C,D,E);CHKERRQ(ierr)
69: #define REG_INTG(A,B,C,D,E) ierr=PetscBagRegisterInt(A,B,C,D,E);CHKERRQ(ierr)
70: #define REG_TRUE(A,B,C,D,E) ierr=PetscBagRegisterBool(A,B,C,D,E);CHKERRQ(ierr)
71: #define REG_STRG(A,B,C,D,E,F) ierr=PetscBagRegisterString(A,B,C,D,E,F);CHKERRQ(ierr)
72: #define REG_ENUM(A,B,C,D,E,F) ierr=PetscBagRegisterEnum(A,B,C,D,E,F);CHKERRQ(ierr)
74: /*-----------------------------------------------------------------------*/
77: int main(int argc,char **argv)
78: /*-----------------------------------------------------------------------*/
79: {
80: DMMG *dmmg; /* multilevel grid structure */
81: AppCtx *user; /* user-defined work context */
82: Parameter *param;
83: int ierr,result = 0;
84: MPI_Comm comm;
85: DM da;
87: PetscInitialize(&argc,&argv,(char *)0,help);
88: comm = PETSC_COMM_WORLD;
90: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
91: Set up the problem parameters.
92: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
93: PetscMalloc(sizeof(AppCtx),&user);
94: PetscBagCreate(comm,sizeof(Parameter),&(user->bag));
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,user->grid.mglevels,user,&dmmg);
104: DMDACreate2d(comm,user->grid.periodic,user->grid.periodic,user->grid.stencil,user->grid.ni,user->grid.nj,PETSC_DECIDE,PETSC_DECIDE,user->grid.dof,user->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:
121: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
122: Free work space.
123: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
124: DMRestoreGlobalVector(da, &(user->Xold));
125: PetscBagDestroy(&user->bag);
126: PetscFree(user);
127: DMDestroy(&da);
128: DMMGDestroy(dmmg);
129: PetscFinalize();
130: return result;
131: }
133: /*---------------------------------------------------------------------*/
136: int SetParams(AppCtx *user)
137: /*---------------------------------------------------------------------*/
138: {
139: PetscBag bag = user->bag;
140: Parameter *p;
141: int ierr, ierr_out=0;
142: PetscBagGetData(bag,(void**)&p);
144: /* give the bag a name */
145: PetscBagSetName(bag,"ex1_params","Parameter bag for ex1.c");
147: /* domain geometry & grid size */
148: REG_INTG(bag,&p->ni,40 ,"ni","Grid points in x-dir");
149: REG_INTG(bag,&p->nj,40 ,"nj","Grid points in y-dir");
150: user->grid.dx = 1.0/((double)(p->ni - 1));
151: user->grid.dz = 1.0/((double)(p->nj - 1));
152: user->grid.ni = p->ni;
153: user->grid.nj = p->nj;
155: /* initial conditions */
156: REG_INTG(bag,&p->flow_type,SHEAR_CELL ,"flow_type","Flow field mode: 0=shear cell, 1=translation");
157: REG_REAL(bag,&p->sigma,0.07 ,"sigma","Standard deviation of the gaussian IC");
158: REG_REAL(bag,&p->xctr,0.5 ,"xctr","x-position of the center of the gaussian IC");
159: REG_REAL(bag,&p->zctr,0.75 ,"zctr","z-position of the center of the gaussian IC");
161: /* time stepping */
162: REG_REAL(bag,&p->t_max,1 ,"t_max","Maximum dimensionless time");
163: REG_REAL(bag,&p->cfl,5 ,"cfl","Courant number");
164: REG_REAL(bag,&p->t_output_interval,0.1 ,"t_output","Dimensionless time interval for output");
165: REG_INTG(bag,&p->N_steps,1000 ,"nsteps","Maximum time-steps");
166: REG_INTG(bag,&p->n,1 ,"nsteps","<DO NOT SET> current time-step");
167: REG_REAL(bag,&p->t,0.0 ,"time","<DO NOT SET> initial time");
168: REG_REAL(bag,&p->dt,p->cfl*PetscMin(user->grid.dx,user->grid.dz),"dt","<DO NOT SET> time-step size");
170: /* output options */
171: REG_TRUE(bag,&p->param_test ,PETSC_FALSE ,"test","Run parameter test only (T/F)");
172: REG_STRG(bag,&p->output_filename,FNAME_LENGTH ,"null","output_file","Name base for output files, set with: -output_file <filename>");
173: 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");
174: p->output_to_file = OptionsHasName("-output_file");
176: user->grid.periodic = DMDA_BOUNDARY_PERIODIC;
177: user->grid.stencil = DMDA_STENCIL_BOX;
178: user->grid.dof = 6;
179: user->grid.stencil_width = 2;
180: user->grid.mglevels = 1;
181: return ierr_out;
182: }
184: /*---------------------------------------------------------------------*/
187: int ReportParams(AppCtx *user)
188: /*---------------------------------------------------------------------*/
189: {
190: Parameter *param;
191: GridInfo grid = user->grid;
192: int ierr, ierr_out=0;
193: PetscBagGetData(user->bag,(void**)¶m);
195: PetscPrintf(PETSC_COMM_WORLD,"---------------MOC test 1----------------\n");
196: PetscPrintf(PETSC_COMM_WORLD,"Prescribed wind, method of\n");
197: PetscPrintf(PETSC_COMM_WORLD,"characteristics advection, explicit time-\n");
198: PetscPrintf(PETSC_COMM_WORLD,"stepping.\n\n");
199: if (param->flow_type == 0) {
200: PetscPrintf(PETSC_COMM_WORLD,"Flow_type: %d (shear cell).\n\n",param->flow_type);
201: }
202: if (param->flow_type == 1) {
203: PetscPrintf(PETSC_COMM_WORLD,"Flow_type: %d (rigid body rotation).\n\n",param->flow_type);
204: }
205: PetscPrintf(PETSC_COMM_WORLD," [ni,nj] = %d, %d [dx,dz] = %5.4g, %5.4g\n",grid.ni,grid.nj,grid.dx,grid.dz);
206: PetscPrintf(PETSC_COMM_WORLD," t_max = %g, cfl = %g, dt = %5.4g,",param->t_max,param->cfl,param->dt);
207: PetscPrintf(PETSC_COMM_WORLD," t_output = %g\n",param->t_output_interval);
208: if (param->output_to_file) {
209: PetscPrintf(PETSC_COMM_WORLD,"Output File: Binary file \"%s\"\n",param->output_filename);
210: }
211: if (!param->output_to_file)
212: PetscPrintf(PETSC_COMM_WORLD,"Output File: NO OUTPUT!\n");
213: PetscPrintf(PETSC_COMM_WORLD,"----------------------------------------\n");
214: if (param->param_test) PetscEnd();
215: return ierr_out;
216: }
218: /* ------------------------------------------------------------------- */
221: int Initialize(DMMG *dmmg)
222: /* ------------------------------------------------------------------- */
223: {
224: AppCtx *user = (AppCtx*)dmmg[0]->user;
225: Parameter *param;
226: DM da;
227: PetscReal PI = 3.14159265358979323846;
228: PetscReal sigma,xc,zc;
229: PetscReal dx=user->grid.dx,dz=user->grid.dz;
230: int i,j,ierr,is,js,im,jm;
231: Field **x;
232: PetscBagGetData(user->bag,(void**)¶m);
233: sigma=param->sigma; xc=param->xctr; zc=param->zctr;
235: /* Get the DM and grid */
236: da = (dmmg[0]->dm);
237: DMDAGetCorners(da,&is,&js,PETSC_NULL,&im,&jm,PETSC_NULL);
238: DMDAVecGetArray(da,user->Xold,(void**)&x);
240: for (j=js; j<js+jm; j++) {
241: for (i=is; i<is+im; i++) {
242: if (param->flow_type == SHEAR_CELL) {
243: x[j][i].u = -sin(PI*i*dx)*cos(PI*j*dz)/dx;
244: x[j][i].w = sin(PI*j*dz)*cos(PI*i*dx)/dz;
245: } else {
246: x[j][i].u = 0.0;
247: x[j][i].w = -1.0/dz;
248: }
249: x[j][i].c1 = 100*exp(-0.5*((i*dx-xc)*(i*dx-xc)+(j*dz-zc)*(j*dz-zc))/sigma/sigma);
250: x[j][i].c2 = x[j][i].c3 = x[j][i].c4 = x[j][i].c1;
251: }
252: }
253:
254: /* restore the grid to it's vector */
255: DMDAVecRestoreArray(da,user->Xold,(void**)&x);
256: VecCopy(user->Xold, DMMGGetx(dmmg));
257: return 0;
258: }
260: /* ------------------------------------------------------------------- */
263: int DoSolve(DMMG *dmmg)
264: /* ------------------------------------------------------------------- */
265: {
266: AppCtx *user = (AppCtx*)dmmg[0]->user;
267: Parameter *param;
268: PetscReal t_output = 0.0;
269: int ierr, i, n_plot = 0, Ncomponents, components[5];
270: DM da = DMMGGetDM(dmmg);
271: Vec Xstar;
272: Characteristic c;
273: PetscBagGetData(user->bag,(void**)¶m);
275: DMGetGlobalVector(da, &Xstar);
277: /*------------ BEGIN CHARACTERISTIC SETUP ---------------*/
278: CharacteristicCreate(PETSC_COMM_WORLD, &c);
279: /* set up the velocity interpolation system */
280: Ncomponents = 2; for (i=0;i<Ncomponents;i++) {components[i] = i;}
281: CharacteristicSetVelocityInterpolationLocal(c, da, DMMGGetx(dmmg), user->Xold, Ncomponents, components, InterpVelocity2D, user);
282: /* set up the fields interpolation system */
283: Ncomponents = 4; for (i=0;i<Ncomponents;i++) {components[i] = i+2;}
284: CharacteristicSetFieldInterpolationLocal(c, da, user->Xold, Ncomponents, components, InterpFields2D, user);
285: /*------------ END CHARACTERISTIC SETUP ----------------*/
287: /* output initial data */
288: PetscPrintf(PETSC_COMM_WORLD," Initialization, Time: %5.4g\n", param->t);
289: DoOutput(dmmg,n_plot);
290: t_output += param->t_output_interval; n_plot++;
292: /* timestep loop */
293: for (param->t=param->dt; param->t<=param->t_max; param->t+=param->dt) {
294: if (param->n > param->N_steps) {
295: PetscPrintf(PETSC_COMM_WORLD,"EXCEEDED MAX NUMBER OF TIMESTEPS! EXITING SOLVE!\n");
296: return 0;
297: }
299: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
300: Solve at time t & copy solution into solution vector.
301: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
302: /* Copy in the velocities to Xstar */
303: VecCopy(DMMGGetx(dmmg), Xstar);
304: /* Put \phi_* into Xstar */
305: CharacteristicSolve(c, param->dt, Xstar);
306: /* Copy the advected field into the solution \phi_t = \phi_* */
307: VecCopy(Xstar, DMMGGetx(dmmg));
309: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
310: Copy new solution to old solution in prep for the next timestep.
311: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
312: VecCopy(DMMGGetx(dmmg), user->Xold);
314: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
315: Timestep complete, report and update counter.
316: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
317: PetscPrintf(PETSC_COMM_WORLD," Step: %d, Time: %5.4g\n", param->n, param->t);
318: param->n++;
320: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
321: Make output.
322: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
323: if (param->t >= t_output) {
324: DoOutput(dmmg,n_plot);
325: t_output += param->t_output_interval; n_plot++;
326: }
327: }
328: DMRestoreGlobalVector(da, &Xstar);
329: CharacteristicDestroy(&c);
330: return 0;
331: }
333: /*---------------------------------------------------------------------*/
336: /* uses analytic velocity fields */
337: PetscErrorCode InterpVelocity2D(void *f, PetscReal ij_real[], PetscInt numComp,
338: PetscInt components[], PetscReal velocity[],
339: void *ctx)
340: /*---------------------------------------------------------------------*/
341: {
342: AppCtx *user = (AppCtx *) ctx;
343: Parameter *param;
344: PetscReal dx=user->grid.dx, dz=user->grid.dz;
345: PetscReal PI = 3.14159265358979323846;
346: int ierr;
347: PetscBagGetData(user->bag,(void**)¶m);
349: if (param->flow_type == SHEAR_CELL) {
350: velocity[0] = -sin(PI*ij_real[0]*dx)*cos(PI*ij_real[1]*dz)/dx;
351: velocity[1] = sin(PI*ij_real[1]*dz)*cos(PI*ij_real[0]*dx)/dz;
352: } else {
353: velocity[0] = 0.;
354: velocity[1] = -1./dz;
355: }
356: return 0;
357: }
359: /*---------------------------------------------------------------------*/
362: /* uses bicubic interpolation */
363: PetscErrorCode InterpFields2D(void *f, PetscReal ij_real[], PetscInt numComp,
364: PetscInt components[], PetscReal field[],
365: void *ctx)
366: /*---------------------------------------------------------------------*/
367: {
368: AppCtx *user = (AppCtx*)ctx;
369: FieldArr **x = (FieldArr**)f;
370: int c, ni=user->grid.ni, nj=user->grid.nj;
371: PetscReal ir=ij_real[0], jr=ij_real[1];
373: for (c=0;c<numComp;c++) {
374: /* boundary condition: set to zero if characteristic begins outside the domain */
375: if ( ir < 0 || ir > ni-1 || jr < 0 || jr> nj-1 ) {
376: field[c] = 0.0;
377: } else {
378: field[c] = BiCubicInterp(x, ir, jr, components[c]);
379: }
380: }
381: return 0;
382: }
384: /*---------------------------------------------------------------------*/
387: PetscReal BiCubicInterp(FieldArr **x, PetscReal ir, PetscReal jr, int c)
388: /*---------------------------------------------------------------------*/
389: {
390: int im, jm, imm,jmm,ip,jp,ipp,jpp;
391: PetscReal il, jl, row1, row2, row3, row4;
392: im = (int)floor(ir); jm = (int)floor(jr);
393: il = ir - im + 1.0; jl = jr - jm + 1.0;
394: imm = im-1; ip = im+1; ipp = im+2;
395: jmm = jm-1; jp = jm+1; jpp = jm+2;
396: row1 = CubicInterp(il,x[jmm][imm][c],x[jmm][im][c],x[jmm][ip][c],x[jmm][ipp][c]);
397: row2 = CubicInterp(il,x[jm] [imm][c],x[jm] [im][c],x[jm] [ip][c],x[jm] [ipp][c]);
398: row3 = CubicInterp(il,x[jp] [imm][c],x[jp] [im][c],x[jp] [ip][c],x[jp] [ipp][c]);
399: row4 = CubicInterp(il,x[jpp][imm][c],x[jpp][im][c],x[jpp][ip][c],x[jpp][ipp][c]);
400: return CubicInterp(jl,row1,row2,row3,row4);
401: }
403: /*---------------------------------------------------------------------*/
406: PetscReal CubicInterp(PetscReal x, PetscReal y_1, PetscReal y_2,
407: PetscReal y_3, PetscReal y_4)
408: /*---------------------------------------------------------------------*/
409: {
410: PetscReal sxth=0.16666666666667, retval;
411: 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
412: - 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;
413: return retval;
414: }
416: /*---------------------------------------------------------------------*/
419: int DoOutput(DMMG *dmmg, int n_plot)
420: /*---------------------------------------------------------------------*/
421: {
422: AppCtx *user = (AppCtx*)dmmg[0]->user;
423: Parameter *param;
424: int ierr;
425: char filename[FNAME_LENGTH];
426: PetscViewer viewer;
427: DM da;
428: PetscBagGetData(user->bag,(void**)¶m);
429: da = DMMGGetDM(dmmg);
431: if (param->output_to_file) { /* send output to binary file */
432: /* generate filename for time t */
433: sprintf(filename,"%s_%3.3d",param->output_filename,n_plot);
434: PetscPrintf(PETSC_COMM_WORLD,"Generating output: time t = %g, ",param->t);
435: PetscPrintf(PETSC_COMM_WORLD,"file = \"%s\"\n",filename);
437: /* make output files */
438: PetscViewerBinaryMatlabOpen(PETSC_COMM_WORLD,filename,&viewer);
439: PetscViewerBinaryMatlabOutputBag(viewer,"par",user->bag);
440: DMDASetFieldNames("u","v","c1","c2","c3","c4",da);
441: PetscViewerBinaryMatlabOutputVecDA(viewer,"field",DMMGGetx(dmmg),da);
442: PetscViewerBinaryMatlabDestroy(&viewer);
443: }
444: return 0;
445: }
447: /* ------------------------------------------------------------------- */
450: int DMDASetFieldNames(const char n0[], const char n1[], const char n2[],
451: const char n3[], const char n4[], const char n5[],
452: DM da)
453: /* ------------------------------------------------------------------- */
454: {
455: int ierr;
456: DMDASetFieldName(da,0,n0);
457: DMDASetFieldName(da,1,n1);
458: DMDASetFieldName(da,2,n2);
459: DMDASetFieldName(da,3,n3);
460: DMDASetFieldName(da,4,n4);
461: DMDASetFieldName(da,5,n5);
462: return 0;
463: }
465: /* ------------------------------------------------------------------- */
468: PetscBool OptionsHasName(const char name[])
469: /* ------------------------------------------------------------------- */
470: {
471: PetscBool retval;
472: int ierr;
473: PetscOptionsHasName(PETSC_NULL,name,&retval);
474: return retval;
475: }