Actual source code: ex4.c
1: static char help[] = "\n";
2: /* - - - - - - - - - - - - - - - - - - - - - - - -
3: ex4.c
4: Advection-diffusion equation
6: Du
7: -- - \nabla^2 u = 0
8: dt
10: which we discretize using Crank-Nicholson
12: u^+ - u^* 1 / \
13: -------- - - | (\nabla^2 u)^+ + (\nabla^2 u)^* | = 0
14: \Delta t 2 \ /
16: which gives, collecting terms
18: \Delta t \Delta t
19: u^+ - -------- (\nabla^2 u)^+ = u^* + -------- (\nabla^2 u)^*
20: 2 2
22: where u^- is the solution at the old time, u^+ is
23: the solution at the new time, and u^* is the solution
24: at the old time advected to the new time.
25: - - - - - - - - - - - - - - - - - - - - - - - - */
26: #include <petscsnes.h>
27: #include <petscdmda.h>
28: #include <petscdmmg.h>
29: #include <petscbag.h>
30: #include <petsccharacteristic.h>
32: #define EXAMPLE_NUMBER 1
33: #define SHEAR_CELL 0
34: #define SOLID_BODY 1
35: #define FNAME_LENGTH 60
37: typedef struct field_s {
38: PetscReal phi;
39: } Field;
41: typedef struct parameter_s {
42: int ni, nj,pi,pj; /* number of grid points, number of procs */
43: PassiveScalar amp,sigma,xctr,zctr,L1,L2,LINF; /* parameters for gaussian Initial Condition */
44: PassiveScalar Pe, theta, ct, st, diffScale; /* parameters for velocity field and diffusion length */
45: int flow_type, sl_event;
46: PetscBool param_test, output_to_file;
47: char output_filename[FNAME_LENGTH];
48: /* timestep stuff */
49: PassiveScalar t; /* the time */
50: int n; /* the time step */
51: PassiveScalar dtAdvection, dtDiffusion; /* minimal advection and diffusion time steps */
52: PassiveScalar t_max, dt, cfl, t_output_interval;
53: int N_steps;
54: } Parameter;
56: typedef struct gridinfo_s {
57: DMDABoundaryType periodic;
58: DMDAStencilType stencil;
59: int ni,nj,dof,stencil_width,mglevels;
60: PassiveScalar dx,dz;
61: } GridInfo;
63: typedef struct appctx_s {
64: DMMG *dmmg;
65: Vec Xold;
66: PetscBag bag;
67: GridInfo *grid;
68: } AppCtx;
70: /* Main routines */
71: int SetParams (AppCtx*);
72: int ReportParams (AppCtx*);
73: int Initialize (DMMG*);
74: int DoSolve (DMMG*);
75: int DoOutput (DMMG*, int);
76: PetscReal BiCubicInterp (Field**, PetscReal, PetscReal);
77: PetscReal CubicInterp (PetscReal, PetscReal, PetscReal, PetscReal, PetscReal);
78: PetscBool OptionsHasName(const char*);
80: /* characteristic call-backs (static interface) */
81: PetscErrorCode InterpVelocity2D(void*, PetscReal[], PetscInt, PetscInt[], PetscReal[], void*);
82: PetscErrorCode InterpFields2D (void*, PetscReal[], PetscInt, PetscInt[], PetscReal[], void*);
84: PetscErrorCode FormOldTimeFunctionLocal(DMDALocalInfo *, PetscScalar **, PetscScalar **, AppCtx *);
85: PetscErrorCode FormNewTimeFunctionLocal(DMDALocalInfo *, PetscScalar **, PetscScalar **, AppCtx *);
87: /* a few macros for convenience */
88: #define REG_REAL(A,B,C,D,E) ierr=PetscBagRegisterReal(A,B,C,D,E);CHKERRQ(ierr)
89: #define REG_INTG(A,B,C,D,E) ierr=PetscBagRegisterInt(A,B,C,D,E);CHKERRQ(ierr)
90: #define REG_TRUE(A,B,C,D,E) ierr=PetscBagRegisterBool(A,B,C,D,E);CHKERRQ(ierr)
91: #define REG_STRG(A,B,C,D,E,F) ierr=PetscBagRegisterString(A,B,C,D,E,F);CHKERRQ(ierr)
92: #define REG_ENUM(A,B,C,D,E,F) ierr=PetscBagRegisterEnum(A,B,C,D,E,F);CHKERRQ(ierr)
94: /*-----------------------------------------------------------------------*/
97: int main(int argc,char **argv)
98: /*-----------------------------------------------------------------------*/
99: {
100: AppCtx *user; /* user-defined work context */
101: Parameter *param;
102: DM da;
103: GridInfo grid;
104: MPI_Comm comm;
105: PetscErrorCode ierr;
107: PetscInitialize(&argc,&argv,(char *)0,help);
108: comm = PETSC_COMM_WORLD;
110: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
111: Set up the problem parameters.
112: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
113: PetscMalloc(sizeof(AppCtx),&user);
114: PetscBagCreate(comm,sizeof(Parameter),&(user->bag));
115: user->grid = &grid;
116: SetParams(user);
117: ReportParams(user);
118: PetscBagGetData(user->bag,(void**)¶m);
120: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
121: Create distributed array multigrid object (DMMG) to manage parallel grid and vectors
122: for principal unknowns (x) and governing residuals (f)
123: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
124: DMMGCreate(comm,grid.mglevels,user,&user->dmmg);
125: DMDACreate2d(comm,grid.periodic,grid.periodic,grid.stencil,grid.ni,grid.nj,PETSC_DECIDE,PETSC_DECIDE,grid.dof,grid.stencil_width,0,0,&da);
126: DMMGSetDM(user->dmmg,(DM) da);
127: DMDestroy(&da);
128: DMMGSetSNESLocal(user->dmmg,FormNewTimeFunctionLocal,PETSC_NULL,PETSC_NULL,PETSC_NULL);
129: DMMGSetFromOptions(user->dmmg);
130: DMDAGetInfo(DMMGGetDM(user->dmmg),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);
131: REG_INTG(user->bag,¶m->pi,param->pi ,"procs_x","<DO NOT SET> Processors in the x-direction");
132: REG_INTG(user->bag,¶m->pj,param->pj ,"procs_y","<DO NOT SET> Processors in the y-direction");
134: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
135: Create user context, set problem data, create vector data structures.
136: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
137: DMGetGlobalVector(DMMGGetDM(user->dmmg), &(user->Xold));
139: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
140: Initialize and solve the nonlinear system
141: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
142: Initialize(user->dmmg);
143: DoSolve(user->dmmg);
144:
145: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
146: Free work space.
147: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
148: DMRestoreGlobalVector(DMMGGetDM(user->dmmg), &(user->Xold));
149: PetscBagDestroy(&user->bag);
150: DMMGDestroy(user->dmmg);
151: PetscFree(user);
152: PetscFinalize();
153: return 0;
154: }
156: /*---------------------------------------------------------------------*/
159: int SetParams(AppCtx *user)
160: /*---------------------------------------------------------------------*/
161: {
162: PetscBag bag = user->bag;
163: Parameter *p;
164: GridInfo *grid = user->grid;
165: int ierr, ierr_out=0;
166: PetscReal PI = 3.14159265358979323846;
168: PetscBagGetData(bag,(void**)&p);
170: /* domain geometry & grid size */
171: REG_INTG(bag,&p->ni,40 ,"ni","Grid points in x-dir");
172: REG_INTG(bag,&p->nj,40 ,"nj","Grid points in y-dir");
173: grid->ni = p->ni;
174: grid->nj = p->nj;
175: grid->dx = 1.0/((double)(grid->ni - 1));
176: grid->dz = 1.0/((double)(grid->nj - 1));
178: /* initial conditions */
179: REG_INTG(bag,&p->flow_type,SOLID_BODY ,"flow_type","Flow field mode: 0=shear cell, 1=translation");
180: REG_REAL(bag,&p->amp,1 ,"amp","Amplitude of the gaussian IC");
181: REG_REAL(bag,&p->sigma,0.07 ,"sigma","Standard deviation of the gaussian IC");
182: REG_REAL(bag,&p->xctr,0.5 ,"xctr","x-position of the center of the gaussian IC");
183: REG_REAL(bag,&p->zctr,0.5 ,"zctr","z-position of the center of the gaussian IC");
185: /* Velocity field */
186: REG_REAL(bag,&p->Pe,10 ,"Pe","Peclet number");
187: REG_REAL(bag,&p->theta,-90 ,"theta","velocity angle (degrees)");
188: REG_REAL(bag,&p->ct,cos(p->theta*PI/180),"cosTheta","<DO NOT SET> cosine velocity angle");
189: REG_REAL(bag,&p->st,sin(p->theta*PI/180),"sinTheta","<DO NOT SET> sine velocity angle");
190:
191: /* diffusion LengthScale for time stepping */
192: REG_REAL(bag,&p->diffScale,2, "diffScale","diffusion length scale (number of grid points for stable diffusion)");
194: /* time stepping */
195: REG_REAL(bag,&p->t_max,1 ,"t_max","Maximum dimensionless time");
196: REG_REAL(bag,&p->cfl,5 ,"cfl","Courant number");
197: REG_REAL(bag,&p->t_output_interval,0.1 ,"t_output","Dimensionless time interval for output");
198: REG_INTG(bag,&p->N_steps,1000 ,"nsteps","Maximum time-steps");
199: REG_INTG(bag,&p->n,1 ,"nsteps","<DO NOT SET> current time-step");
200: REG_REAL(bag,&p->t,0.0 ,"time","<DO NOT SET> initial time");
201: REG_REAL(bag,&p->dtAdvection,1./p->Pe/(sqrt(pow(p->ct/grid->dx,2)+pow(p->st/grid->dz,2))),"dtAdvection","<DO NOT SET> CFL limited advection time step");
202: REG_REAL(bag,&p->dtDiffusion,1./(1./grid->dx/grid->dx+1./grid->dz/grid->dz),"dtDiffusion","<DO NOT SET> grid-space limited diffusion time step");
203: REG_REAL(bag,&p->dt,PetscMin(p->cfl*p->dtAdvection,pow(p->diffScale,2)*p->dtDiffusion),"dt","<DO NOT SET> time-step size");
205: /* output options */
206: REG_TRUE(bag,&p->param_test ,PETSC_FALSE ,"test","Run parameter test only (T/F)");
207: REG_STRG(bag,&p->output_filename,FNAME_LENGTH ,"null","output_file","Name base for output files, set with: -output_file <filename>");
208: 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");
209: p->output_to_file = OptionsHasName("-output_file");
211: grid->ni = p->ni;
212: grid->nj = p->nj;
213: grid->periodic = DMDA_BOUNDARY_PERIODIC;
214: grid->stencil = DMDA_STENCIL_BOX;
215: grid->dof = 1;
216: grid->stencil_width = 2;
217: grid->mglevels = 1;
218: return ierr_out;
219: }
221: /*---------------------------------------------------------------------*/
224: int ReportParams(AppCtx *user)
225: /*---------------------------------------------------------------------*/
226: {
227: Parameter *param;
228: GridInfo *grid = user->grid;
229: int ierr, ierr_out=0;
230: PetscBagGetData(user->bag,(void**)¶m);
232: PetscPrintf(PETSC_COMM_WORLD,"---------------MOC test 1----------------\n");
233: PetscPrintf(PETSC_COMM_WORLD,"Prescribed wind, method of\n");
234: PetscPrintf(PETSC_COMM_WORLD,"characteristics advection, explicit time-\n");
235: PetscPrintf(PETSC_COMM_WORLD,"stepping.\n\n");
236: if (param->flow_type == 0) {
237: PetscPrintf(PETSC_COMM_WORLD,"Flow_type: %d (shear cell).\n\n",param->flow_type);
238: }
239: if (param->flow_type == 1) {
240: PetscPrintf(PETSC_COMM_WORLD,"Flow_type: %d (translation).\n\n",param->flow_type);
241: }
242: PetscPrintf(PETSC_COMM_WORLD," [ni,nj] = %d, %d [dx,dz] = %5.4g, %5.4g\n",grid->ni,grid->nj,grid->dx,grid->dz);
243: PetscPrintf(PETSC_COMM_WORLD," t_max = %g, cfl = %g, dt = %5.4g,",param->t_max,param->cfl,param->dt);
244: PetscPrintf(PETSC_COMM_WORLD," t_output = %g\n",param->t_output_interval);
245: PetscPrintf(PETSC_COMM_WORLD," dt_advection= %g, dt_diffusion= %g\n",param->dtAdvection,param->dtDiffusion);
246: if (param->output_to_file) {
247: PetscPrintf(PETSC_COMM_WORLD,"Output File: Binary file \"%s\"\n",param->output_filename);
248: }
249: if (!param->output_to_file)
250: PetscPrintf(PETSC_COMM_WORLD,"Output File: NO OUTPUT!\n");
251: PetscPrintf(PETSC_COMM_WORLD,"----------------------------------------\n");
252: if (param->param_test) PetscEnd();
253: return ierr_out;
254: }
256: /* ------------------------------------------------------------------- */
259: int Initialize(DMMG *dmmg)
260: /* ------------------------------------------------------------------- */
261: {
262: AppCtx *user = (AppCtx*)dmmg[0]->user;
263: Parameter *param;
264: DM da;
265: PetscReal amp, sigma, xc, zc ;
266: PetscReal dx=user->grid->dx,dz=user->grid->dz;
267: int i,j,ierr,is,js,im,jm;
268: Field **x;
269: PetscBagGetData(user->bag,(void**)¶m);
270:
271: amp = param->amp;
272: sigma = param->sigma;
273: xc = param->xctr; zc = param->zctr;
275: /* Get the DM and grid */
276: da = DMMGGetDM(dmmg);
277: DMDAGetCorners(da,&is,&js,PETSC_NULL,&im,&jm,PETSC_NULL);
278: DMDAVecGetArray(da,user->Xold,(void**)&x);
280: for (j=js; j<js+jm; j++) {
281: for (i=is; i<is+im; i++) {
282: x[j][i].phi = param->amp*exp(-0.5*((i*dx-xc)*(i*dx-xc)+(j*dz-zc)*(j*dz-zc))/sigma/sigma);
283: }
284: }
285:
286: /* restore the grid to it's vector */
287: DMDAVecRestoreArray(da,user->Xold,(void**)&x);
288: VecCopy(user->Xold, DMMGGetx(dmmg));
290: return 0;
291: }
293: /* ------------------------------------------------------------------- */
296: int DoSolve(DMMG *dmmg)
297: /* ------------------------------------------------------------------- */
298: {
299: AppCtx *user = (AppCtx*)dmmg[0]->user;
300: Parameter *param;
301: PetscReal t_output = 0.0;
302: int ierr, n_plot = 0, Ncomponents, components[3];
303: DM da = DMMGGetDM(dmmg);
304: Vec Xstar;
305: Characteristic c;
306: PetscBagGetData(user->bag,(void**)¶m);
308: DMGetGlobalVector(da, &Xstar);
310: /*------------ BEGIN CHARACTERISTIC SETUP ---------------*/
311: CharacteristicCreate(PETSC_COMM_WORLD, &c);
312: /* set up the velocity interpolation system */
313: Ncomponents = 2; components[0] = 0; components[1] = 0;
314: CharacteristicSetVelocityInterpolationLocal(c, da, DMMGGetx(dmmg), user->Xold, Ncomponents, components, InterpVelocity2D, user);
315: /* set up the fields interpolation system */
316: Ncomponents = 1; components[0] = 0;
317: CharacteristicSetFieldInterpolationLocal(c, da, user->Xold, Ncomponents, components, InterpFields2D, user);
318: /*------------ END CHARACTERISTIC SETUP ----------------*/
320: /* output initial data */
321: PetscPrintf(PETSC_COMM_WORLD," Initialization, Time: %5.4g\n", param->t);
322: DoOutput(dmmg,n_plot);
323: t_output += param->t_output_interval; n_plot++;
325: /* timestep loop */
326: for (param->t=param->dt; param->t<=param->t_max; param->t+=param->dt) {
327: if (param->n > param->N_steps) {
328: PetscPrintf(PETSC_COMM_WORLD,"EXCEEDED MAX NUMBER OF TIMESTEPS! EXITING SOLVE!\n");
329: return 0;
330: }
332: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
333: Solve at time t & copy solution into solution vector.
334: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
335: /* Evaluate operator (I + \Delta t/2 L) u^- = X^- */
336: DMDAFormFunctionLocal(da, (DMDALocalFunction1) FormOldTimeFunctionLocal, DMMGGetx(dmmg), user->Xold, user);
337: /* Advect Xold into Xstar */
338: CharacteristicSolve(c, param->dt, Xstar);
339: /* Xstar -> Xold */
340: VecCopy(Xstar, user->Xold);
341: /* Solve u^+ = (I - \Delta t/2 L)^-1 Xstar which could be F(u^+) = Xstar */
342: DMMGSolve(dmmg);
344: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
345: report step and update counter.
346: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
347: PetscPrintf(PETSC_COMM_WORLD," Step: %d, Time: %5.4g\n", param->n, param->t);
348: param->n++;
350: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
351: Output variables.
352: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
353: if (param->t >= t_output) {
354: DoOutput(dmmg,n_plot);
355: t_output += param->t_output_interval; n_plot++;
356: }
357: }
358: DMRestoreGlobalVector(da, &Xstar);
359: CharacteristicDestroy(&c);
360: return 0;
361: }
363: /*---------------------------------------------------------------------*/
366: /* linear interpolation, ir: [0, ni], jr: [0, nj] */
367: PetscErrorCode InterpVelocity2D(void *f, PetscReal ij_real[], PetscInt numComp,
368: PetscInt components[], PetscReal velocity[],
369: void *ctx)
370: /*---------------------------------------------------------------------*/
371: {
372: AppCtx *user = (AppCtx *) ctx;
373: Parameter *param;
374: PetscReal dx=user->grid->dx, dz=user->grid->dz;
375: PetscReal PI = 3.14159265358979323846;
376: int ierr;
377: PetscBagGetData(user->bag,(void**)¶m);
379: /* remember: must be coordinate velocities not true velocities */
380: if (param->flow_type == SHEAR_CELL) {
381: velocity[0] = -sin(PI*ij_real[0]*dx)*cos(PI*ij_real[1]*dz)/dx;
382: velocity[1] = sin(PI*ij_real[1]*dz)*cos(PI*ij_real[0]*dx)/dz;
383: } else {
384: velocity[0] = param->Pe*param->ct/dx;
385: velocity[1] = param->Pe*param->st/dz;
386: }
387: return 0;
388: }
390: /*---------------------------------------------------------------------*/
393: PetscErrorCode InterpFields2D(void *f, PetscReal ij_real[], PetscInt numComp,
394: PetscInt components[], PetscReal field[],
395: void *ctx)
396: /*---------------------------------------------------------------------*/
397: {
398: AppCtx *user = (AppCtx*)ctx;
399: Field **x = (Field**)f;
400: int ni=user->grid->ni, nj=user->grid->nj;
401: int ierr;
402: PetscReal ir=ij_real[0], jr=ij_real[1];
404: /* map back to periodic domain if out of bounds */
405: if ( ir < 0 || ir > ni-1 || jr < 0 || jr> nj-1 ) {
406: DMDAMapCoordsToPeriodicDomain(DMMGGetDM(user->dmmg), &ir, &jr);
407: }
408: field[0] = BiCubicInterp(x, ir, jr);
409: return 0;
410: }
412: /*---------------------------------------------------------------------*/
415: PetscReal BiCubicInterp(Field **x, PetscReal ir, PetscReal jr)
416: /*---------------------------------------------------------------------*/
417: {
418: int im, jm, imm,jmm,ip,jp,ipp,jpp;
419: PetscReal il, jl, row1, row2, row3, row4;
420: im = (int)floor(ir); jm = (int)floor(jr);
421: il = ir - im + 1.0; jl = jr - jm + 1.0;
422: imm = im-1; ip = im+1; ipp = im+2;
423: jmm = jm-1; jp = jm+1; jpp = jm+2;
424: row1 = CubicInterp(il,x[jmm][imm].phi,x[jmm][im].phi,x[jmm][ip].phi,x[jmm][ipp].phi);
425: row2 = CubicInterp(il,x[jm] [imm].phi,x[jm] [im].phi,x[jm] [ip].phi,x[jm] [ipp].phi);
426: row3 = CubicInterp(il,x[jp] [imm].phi,x[jp] [im].phi,x[jp] [ip].phi,x[jp] [ipp].phi);
427: row4 = CubicInterp(il,x[jpp][imm].phi,x[jpp][im].phi,x[jpp][ip].phi,x[jpp][ipp].phi);
428: return CubicInterp(jl,row1,row2,row3,row4);
429: }
431: /*---------------------------------------------------------------------*/
434: PetscReal CubicInterp(PetscReal x, PetscReal y_1, PetscReal y_2,
435: PetscReal y_3, PetscReal y_4)
436: /*---------------------------------------------------------------------*/
437: {
438: PetscReal sxth=0.16666666666667, retval;
439: 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
440: - 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;
441: return retval;
442: }
444: /*---------------------------------------------------------------------*/
447: int DoOutput(DMMG *dmmg, int n_plot)
448: /*---------------------------------------------------------------------*/
449: {
450: AppCtx *user = (AppCtx*)dmmg[0]->user;
451: Parameter *param;
452: int ierr;
453: char filename[FNAME_LENGTH];
454: PetscViewer viewer;
455: DM da;
456: PetscBagGetData(user->bag,(void**)¶m);
457: da = DMMGGetDM(dmmg);
459: if (param->output_to_file) { /* send output to binary file */
460: /* generate filename for time t */
461: sprintf(filename,"%s_%3.3d",param->output_filename,n_plot);
462: PetscPrintf(PETSC_COMM_WORLD,"Generating output: time t = %g, ",param->t);
463: PetscPrintf(PETSC_COMM_WORLD,"file = \"%s\"\n",filename);
465: /* make output files */
466: PetscViewerBinaryMatlabOpen(PETSC_COMM_WORLD,filename,&viewer);
467: PetscViewerBinaryMatlabOutputBag(viewer,"par",user->bag);
468: DMDASetFieldName(da,0,"phi");
469: PetscViewerBinaryMatlabOutputVecDA(viewer,"field",DMMGGetx(dmmg),da);
470: PetscViewerBinaryMatlabDestroy(&viewer);
471: }
472: return 0;
473: }
475: /* ------------------------------------------------------------------- */
478: PetscBool OptionsHasName(const char name[])
479: /* ------------------------------------------------------------------- */
480: {
481: PetscBool retval;
482: int ierr;
483: PetscOptionsHasName(PETSC_NULL,name,&retval);
484: return retval;
485: }
489: /*
490: FormNewTimeFunctionLocal - Evaluates f = (I - \Delta t/2 L) x - f(u^-)^*.
492: Note: We get f(u^-)^* from Xold in the user context.
494: Process adiC(36): FormNewTimeFunctionLocal
495: */
496: PetscErrorCode FormNewTimeFunctionLocal(DMDALocalInfo *info, PetscScalar **x, PetscScalar **f, AppCtx *user)
497: {
498: DM da = DMMGGetDM(user->dmmg);
499: PetscScalar **fold;
500: Parameter *param;
501: PetscScalar u,uxx,uyy;
502: PetscReal hx,hy,hxdhy,hydhx;
503: PetscInt i,j;
507: hx = 1.0/(PetscReal)(info->mx-1);
508: hy = 1.0/(PetscReal)(info->my-1);
509: hxdhy = hx/hy;
510: hydhx = hy/hx;
512: PetscBagGetData(user->bag, (void**) ¶m);
513: DMDAVecGetArray(da, user->Xold, &fold);
514: for (j=info->ys; j<info->ys+info->ym; j++) {
515: for (i=info->xs; i<info->xs+info->xm; i++) {
516: u = x[j][i];
517: uxx = (2.0*u - x[j][i-1] - x[j][i+1])*hydhx;
518: uyy = (2.0*u - x[j-1][i] - x[j+1][i])*hxdhy;
519: f[j][i] = u*hx*hy + param->dt*0.5*(uxx + uyy) - fold[j][i];
520: }
521: }
522: DMDAVecRestoreArray(da, user->Xold, &fold);
524: PetscLogFlops(13.0*info->ym*info->xm);
525: return(0);
526: }
530: /*
531: FormOldTimeFunctionLocal - Evaluates f = (I + \Delta t/2 L) x.
533: Process adiC(36): FormOldTimeFunctionLocal
534: */
535: PetscErrorCode FormOldTimeFunctionLocal(DMDALocalInfo *info, PetscScalar **x, PetscScalar **f, AppCtx *user)
536: {
537: Parameter *param;
538: PetscScalar u,uxx,uyy;
539: PetscReal hx,hy,hxdhy,hydhx;
540: PetscInt i,j;
544: hx = 1.0/(PetscReal)(info->mx-1);
545: hy = 1.0/(PetscReal)(info->my-1);
546: hxdhy = hx/hy;
547: hydhx = hy/hx;
549: PetscBagGetData(user->bag, (void**) ¶m);
550: for (j=info->ys; j<info->ys+info->ym; j++) {
551: for (i=info->xs; i<info->xs+info->xm; i++) {
552: u = x[j][i];
553: uxx = (2.0*u - x[j][i-1] - x[j][i+1])*hydhx;
554: uyy = (2.0*u - x[j-1][i] - x[j+1][i])*hxdhy;
555: f[j][i] = u*hx*hy - param->dt*0.5*(uxx + uyy);
556: }
557: }
559: PetscLogFlops(12.0*info->ym*info->xm);
560: return(0);
561: }