Actual source code: ex32.c
2: static char help[] = "Model multi-physics solver. Modified from ex19.c \n\\n\
3: The 2D driven cavity problem is solved in a velocity-vorticity formulation.\n\
4: The flow can be driven with the lid or with bouyancy or both:\n\
5: -lidvelocity <lid>, where <lid> = dimensionless velocity of lid\n\
6: -grashof <gr>, where <gr> = dimensionless temperature gradent\n\
7: -prandtl <pr>, where <pr> = dimensionless thermal/momentum diffusity ratio\n\
8: -contours : draw contour plots of solution\n\n";
10: /*T
11: Concepts: SNES^solving a system of nonlinear equations (parallel multicomponent example);
12: Concepts: DMDA^using distributed arrays;
13: Concepts: multicomponent
14: Processors: n
15: T*/
17: /* ------------------------------------------------------------------------
19: We thank David E. Keyes for contributing the driven cavity discretization
20: within this example code.
22: This problem is modeled by the partial differential equation system
23:
24: - Lap(U) - Grad_y(Omega) = 0
25: - Lap(V) + Grad_x(Omega) = 0
26: - Lap(Omega) + Div([U*Omega,V*Omega]) - GR*Grad_x(T) = 0
27: - Lap(T) + PR*Div([U*T,V*T]) = 0
29: in the unit square, which is uniformly discretized in each of x and
30: y in this simple encoding.
32: No-slip, rigid-wall Dirichlet conditions are used for [U,V].
33: Dirichlet conditions are used for Omega, based on the definition of
34: vorticity: Omega = - Grad_y(U) + Grad_x(V), where along each
35: constant coordinate boundary, the tangential derivative is zero.
36: Dirichlet conditions are used for T on the left and right walls,
37: and insulation homogeneous Neumann conditions are used for T on
38: the top and bottom walls.
40: A finite difference approximation with the usual 5-point stencil
41: is used to discretize the boundary value problem to obtain a
42: nonlinear system of equations. Upwinding is used for the divergence
43: (convective) terms and central for the gradient (source) terms.
44:
45: The Jacobian can be either
46: * formed via finite differencing using coloring (the default), or
47: * applied matrix-free via the option -snes_mf
48: (for larger grid problems this variant may not converge
49: without a preconditioner due to ill-conditioning).
51: ------------------------------------------------------------------------- */
52: #include <petscsnes.h>
53: #include <petscdmda.h>
54: #include <petscdmcomposite.h>
55: #include <petscdmmg.h>
57: /* User-defined routines and data structure */
58: typedef struct {
59: PetscScalar u,v,omega,temp;
60: } Field;
62: typedef struct {
63: PetscScalar u,v,omega;
64: } Field1;
66: typedef struct {
67: PetscScalar temp;
68: } Field2;
70: typedef struct {
71: PassiveReal lidvelocity,prandtl,grashof; /* physical parameters */
72: PetscBool draw_contours; /* flag - 1 indicates drawing contours */
73: DMMG *dmmg,*dmmg_comp; /* used by MySolutionView() */
74: DMMG *dmmg1,*dmmg2; /* passing objects of sub-physics into the composite physics - used by FormInitialGuessComp() */
75: DM pack;
76: PetscBool COMPOSITE_MODEL;
77: Field **x; /* passing DMMGGetx(dmmg) - final solution of original coupled physics */
78: Field1 **x1; /* passing local ghosted vector array of Physics 1 */
79: Field2 **x2; /* passing local ghosted vector array of Physics 2 */
80: } AppCtx;
96: int main(int argc,char **argv)
97: {
98: DMMG *dmmg,*dmmg1,*dmmg2,*dmmg_comp; /* multilevel grid structure */
99: AppCtx user; /* user-defined work context */
100: PetscInt mx,my,its,dof,nlevels=1;
102: MPI_Comm comm;
103: SNES snes;
104: DM da;
105: PetscBool ViewSolu=PETSC_FALSE,CompSolu=PETSC_TRUE,DoSubPhysics=PETSC_FALSE;
106: Vec solu_true,solu_local;
108: PetscInitialize(&argc,&argv,(char *)0,help);
109: comm = PETSC_COMM_WORLD;
111: PetscOptionsGetInt(PETSC_NULL,"-nlevels",&nlevels,PETSC_NULL);
112: DMMGCreate(comm,nlevels,&user,&dmmg);
114: /*
115: Create distributed array multigrid object (DMMG) to manage parallel grid and vectors
116: for principal unknowns (x) and governing residuals (f)
117: */
118: dof = 4;
119: DMDACreate2d(PETSC_COMM_WORLD,DMDA_BOUNDARY_NONE,DMDA_BOUNDARY_NONE,DMDA_STENCIL_STAR,-4,-4,PETSC_DECIDE,PETSC_DECIDE,dof,1,0,0,&da);
120: DMMGSetDM(dmmg,(DM)da);
121: DMDestroy(&da);
123: DMDAGetInfo(DMMGGetDM(dmmg),0,&mx,&my,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,
124: PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE);
126: /* Problem parameters (velocity of lid, prandtl, and grashof numbers) */
127: user.lidvelocity = 1.0/(mx*my);
128: user.prandtl = 1.0;
129: user.grashof = 1.0;
130: PetscOptionsGetReal(PETSC_NULL,"-lidvelocity",&user.lidvelocity,PETSC_NULL);
131: PetscOptionsGetReal(PETSC_NULL,"-prandtl",&user.prandtl,PETSC_NULL);
132: PetscOptionsGetReal(PETSC_NULL,"-grashof",&user.grashof,PETSC_NULL);
133: PetscOptionsHasName(PETSC_NULL,"-contours",&user.draw_contours);
134: PetscOptionsHasName(PETSC_NULL,"-view_solu",&ViewSolu);
135: PetscOptionsHasName(PETSC_NULL,"-do_subphysics",&DoSubPhysics);
137: DMDASetFieldName(DMMGGetDM(dmmg),0,"x-velocity");
138: DMDASetFieldName(DMMGGetDM(dmmg),1,"y-velocity");
139: DMDASetFieldName(DMMGGetDM(dmmg),2,"Omega");
140: DMDASetFieldName(DMMGGetDM(dmmg),3,"temperature");
142: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
143: Create user context, set problem data, create vector data structures.
144: Also, compute the initial guess.
145: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
147: /* Create nonlinear solver context */
148: DMMGSetSNESLocal(dmmg,FormFunctionLocal,0,ad_FormFunctionLocal,admf_FormFunctionLocal);
149: DMMGSetFromOptions(dmmg);
150: PetscPrintf(comm,"lid velocity = %G, prandtl # = %G, grashof # = %G\n",
151: user.lidvelocity,user.prandtl,user.grashof);
152: DMMGSetInitialGuess(dmmg,FormInitialGuessLocal);
154: /* Solve the nonlinear system */
155: DMMGSolve(dmmg);
157: snes = DMMGGetSNES(dmmg);
158: SNESGetIterationNumber(snes,&its);
159: PetscPrintf(comm,"Origianl Physics: Number of Newton iterations = %D\n\n", its);
160:
161: /* Save the ghosted local solu_true to be used by Physics 1 and Physics 2 */
162: da = DMMGGetDM(dmmg);
163: solu_true = DMMGGetx(dmmg);
164: DMCreateLocalVector(da,&solu_local);
165: DMGlobalToLocalBegin(da,solu_true,INSERT_VALUES,solu_local);
166: DMGlobalToLocalEnd(da,solu_true,INSERT_VALUES,solu_local);
167: DMDAVecGetArray(da,solu_local,(Field **)&user.x);
169: /* Visualize solution */
170: if (user.draw_contours) {
171: VecView(DMMGGetx(dmmg),PETSC_VIEWER_DRAW_WORLD);
172: }
173: if (ViewSolu){
174: user.dmmg = dmmg;
175: MySolutionView(comm,0,&user);
176: }
178: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
179: Setup Physics 1:
180: - Lap(U) - Grad_y(Omega) = 0
181: - Lap(V) + Grad_x(Omega) = 0
182: - Lap(Omega) + Div([U*Omega,V*Omega]) - GR*Grad_x(T) = 0
183: where T is given by the computed x.temp
184: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
185: DMMGCreate(comm,nlevels,&user,&dmmg1);
186: dof = 3;
187: user.COMPOSITE_MODEL = PETSC_FALSE;
188: DMDACreate2d(PETSC_COMM_WORLD,DMDA_BOUNDARY_NONE,DMDA_BOUNDARY_NONE,DMDA_STENCIL_STAR,-4,-4,PETSC_DECIDE,PETSC_DECIDE,dof,1,0,0,&da);
189: DMMGSetDM(dmmg1,(DM)da);
190: DMDASetFieldName(da,0,"x-velocity");
191: DMDASetFieldName(da,1,"y-velocity");
192: DMDASetFieldName(da,2,"Omega");
193: DMDestroy(&da);
195: if (DoSubPhysics){
196: DMMGSetSNESLocal(dmmg1,FormFunctionLocal1,0,ad_FormFunctionLocal,admf_FormFunctionLocal);
197: DMMGSetFromOptions(dmmg1);
198: DMMGSetInitialGuess(dmmg1,FormInitialGuessLocal1);
200: DMMGSolve(dmmg1);
201: snes = DMMGGetSNES(dmmg1);
202: SNESGetIterationNumber(snes,&its);
203: PetscPrintf(comm,"SubPhysics 1, Number of Newton iterations = %D\n\n", its);
204:
205: if (ViewSolu){ /* View individial componets of the solution */
206: user.dmmg1 = dmmg1;
207: MySolutionView(comm,1,&user);
208: }
209: }
211: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
212: Setup Physics 2:
213: - Lap(T) + PR*Div([U*T,V*T]) = 0
214: where U and V are given by the computed x.u and x.v
215: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
216: DMMGCreate(comm,nlevels,&user,&dmmg2);
217: dof = 1;
218: user.COMPOSITE_MODEL = PETSC_FALSE;
219: DMDACreate2d(PETSC_COMM_WORLD,DMDA_BOUNDARY_NONE,DMDA_BOUNDARY_NONE,DMDA_STENCIL_STAR,-4,-4,PETSC_DECIDE,PETSC_DECIDE,dof,1,0,0,&da);
220: DMMGSetDM(dmmg2,(DM)da);
221: DMDASetFieldName(da,0,"temperature");
222: DMDestroy(&da);
224: if (DoSubPhysics){
225: DMMGSetSNESLocal(dmmg2,FormFunctionLocal2,0,ad_FormFunctionLocal,admf_FormFunctionLocal);
226: DMMGSetFromOptions(dmmg2);
227: DMMGSetInitialGuess(dmmg2,FormInitialGuessLocal2);
229: DMMGSolve(dmmg2);
230: snes = DMMGGetSNES(dmmg2);
231: SNESGetIterationNumber(snes,&its);
232: PetscPrintf(comm,"SubPhysics 2, Number of Newton iterations = %D\n\n", its);
234: if (ViewSolu){ /* View individial componets of the solution */
235: user.dmmg2 = dmmg2;
236: MySolutionView(comm,2,&user);
237: }
238: }
239: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
240: Create the DMComposite object to manage the two grids/physics.
241: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
242: DMCompositeCreate(comm,&user.pack);
243: DMCompositeAddDM(user.pack,(DM)DMMGGetDM(dmmg1));
244: DMCompositeAddDM(user.pack,(DM)DMMGGetDM(dmmg2));
246: /* Create the solver object and attach the grid/physics info */
247: DMMGCreate(comm,nlevels,&user,&dmmg_comp);
248: DMMGSetDM(dmmg_comp,(DM)user.pack);
249: CHKMEMQ;
250: DMMGSetISColoringType(dmmg_comp,IS_COLORING_GLOBAL);
252: user.dmmg1 = dmmg1;
253: user.dmmg2 = dmmg2;
254: user.COMPOSITE_MODEL = PETSC_TRUE;
255: DMMGSetInitialGuess(dmmg_comp,FormInitialGuessComp);
256: DMMGSetSNES(dmmg_comp,FormFunctionComp,0);
257: DMMGSetFromOptions(dmmg_comp);
259: /* Solve the nonlinear system */
260: DMMGSolve(dmmg_comp);
261: snes = DMMGGetSNES(dmmg_comp);
262: SNESGetIterationNumber(snes,&its);
263: PetscPrintf(comm,"Composite Physics: Number of Newton iterations = %D\n\n", its);
265: /* Compare the solutions obtained from dmmg and dmmg_comp */
266: if (ViewSolu || CompSolu){
267: /* View the solution of composite physics (phy_num = 3)
268: or Compare solutions (phy_num > 3) */
269: user.dmmg = dmmg;
270: user.dmmg_comp = dmmg_comp;
271: MySolutionView(comm,4,&user);
272: }
274: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
275: Free spaces
276: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
277: DMDAVecRestoreArray(DMMGGetDM(dmmg),solu_local,(Field **)&user.x);
278: VecDestroy(&solu_local);
279: DMDestroy(&user.pack);
280: DMMGDestroy(dmmg);
281: DMMGDestroy(dmmg1);
282: DMMGDestroy(dmmg2);
283: DMMGDestroy(dmmg_comp);
284: PetscFinalize();
285: return 0;
286: }
288: /* ------------------------------------------------------------------- */
292: /*
293: FormInitialGuessLocal - Forms initial approximation for this process
295: Input Parameters:
296: user - user-defined application context
297: X - vector (DMDA local vector)
299: Output Parameter:
300: X - vector with the local values set
301: */
302: PetscErrorCode FormInitialGuessLocal(DMMG dmmg,Vec X)
303: {
304: AppCtx *user = (AppCtx*)dmmg->user;
305: DM da = dmmg->dm;
306: PetscInt i,j,mx,xs,ys,xm,ym;
308: PetscReal grashof,dx;
309: Field **x;
311: grashof = user->grashof;
312: DMDAGetInfo(da,0,&mx,0,0,0,0,0,0,0,0,0,0,0);
313: dx = 1.0/(mx-1);
315: /*
316: Get local grid boundaries (for 2-dimensional DMDA):
317: xs, ys - starting grid indices (no ghost points)
318: xm, ym - widths of local grid (no ghost points)
319: */
320: DMDAGetCorners(da,&xs,&ys,PETSC_NULL,&xm,&ym,PETSC_NULL);
322: /*
323: Get a pointer to vector data.
324: - For default PETSc vectors, VecGetArray() returns a pointer to
325: the data array. Otherwise, the routine is implementation dependent.
326: - You MUST call VecResetArraystoreArray() when you no longer need access to
327: the array.
328: */
329: DMDAVecGetArray(da,X,&x);
331: /*
332: Compute initial guess over the locally owned part of the grid
333: Initial condition is motionless fluid and equilibrium temperature
334: U = V = Omega = 0.0; Temp[x_i, y_j] = x_i;
335: */
336: for (j=ys; j<ys+ym; j++) {
337: for (i=xs; i<xs+xm; i++) {
338: x[j][i].u = 0.0;
339: x[j][i].v = 0.0;
340: x[j][i].omega = 0.0;
341: x[j][i].temp = (grashof>0)*i*dx;
342: }
343: }
345: /* Restore vector */
346: DMDAVecRestoreArray(da,X,&x);
347: return 0;
348: }
352: /* Form initial guess for Physic 1 */
353: PetscErrorCode FormInitialGuessLocal1(DMMG dmmg,Vec X)
354: {
355: DM da = dmmg->dm;
356: PetscInt i,j,mx,xs,ys,xm,ym;
358: Field1 **x;
360: DMDAGetInfo(da,0,&mx,0,0,0,0,0,0,0,0,0,0,0);
361: DMDAGetCorners(da,&xs,&ys,PETSC_NULL,&xm,&ym,PETSC_NULL);
363: DMDAVecGetArray(da,X,&x);
364: for (j=ys; j<ys+ym; j++) {
365: for (i=xs; i<xs+xm; i++) {
366: x[j][i].u = 0.0;
367: x[j][i].v = 0.0;
368: x[j][i].omega = 0.0;
369: }
370: }
371: DMDAVecRestoreArray(da,X,&x);
372: return 0;
373: }
377: /* Form initial guess for Physic 2 */
378: PetscErrorCode FormInitialGuessLocal2(DMMG dmmg,Vec X)
379: {
380: AppCtx *user = (AppCtx*)dmmg->user;
381: DM da = dmmg->dm;
382: PetscInt i,j,mx,xs,ys,xm,ym;
384: PetscReal grashof,dx;
385: Field2 **x;
387: grashof = user->grashof;
388: DMDAGetInfo(da,0,&mx,0,0,0,0,0,0,0,0,0,0,0);
389: dx = 1.0/(mx-1);
391: DMDAGetCorners(da,&xs,&ys,PETSC_NULL,&xm,&ym,PETSC_NULL);
393: DMDAVecGetArray(da,X,&x);
394: for (j=ys; j<ys+ym; j++) {
395: for (i=xs; i<xs+xm; i++) {
396: x[j][i].temp = (grashof>0)*i*dx;
397: }
398: }
399: DMDAVecRestoreArray(da,X,&x);
400: return 0;
401: }
405: /*
406: FormInitialGuessComp -
407: Forms the initial guess for the composite model
408: Unwraps the global solution vector and passes its local pieces into the user functions
409: */
410: PetscErrorCode FormInitialGuessComp(DMMG dmmg,Vec X)
411: {
413: AppCtx *user = (AppCtx*)dmmg->user;
414: DMMG *dmmg1 = user->dmmg1,*dmmg2=user->dmmg2;
415: DM dm = dmmg->dm;
416: Vec X1,X2;
419: /* Access the subvectors in X */
420: DMCompositeGetAccess(dm,X,&X1,&X2);
422: /* Evaluate local user provided function */
423: FormInitialGuessLocal1(*dmmg1,X1);
424: FormInitialGuessLocal2(*dmmg2,X2);
426: DMCompositeRestoreAccess(dm,X,&X1,&X2);
427: return(0);
428: }
432: /*
433: FormFunctionLocal - Function evaluation for this process
435: Input Parameters:
436: info - DMDALocalInfo context
437: x - (DMDA local) vector array including ghost points
438: f - (DMDA local) vector array to be evaluated
439: ptr - user-defined application context
441: Output Parameter:
442: f - array holds local function values
443:
444: f.u = - Lap(U) - Grad_y(Omega)
445: f.v = - Lap(V) + Grad_x(Omega)
446: f.omega = - Lap(Omega) + Div([U*Omega,V*Omega]) - GR*Grad_x(T)
447: f.temp = - Lap(T) + PR*Div([U*T,V*T])
448: */
449: PetscErrorCode FormFunctionLocal(DMDALocalInfo *info,Field **x,Field **f,void *ptr)
450: {
451: AppCtx *user = (AppCtx*)ptr;
453: PetscInt xints,xinte,yints,yinte,i,j;
454: PetscReal hx,hy,dhx,dhy,hxdhy,hydhx;
455: PetscReal grashof,prandtl,lid;
456: PetscScalar u,uxx,uyy,vx,vy,avx,avy,vxp,vxm,vyp,vym;
459: grashof = user->grashof;
460: prandtl = user->prandtl;
461: lid = user->lidvelocity;
463: /*
464: Define mesh intervals ratios for uniform grid.
466: Note: FD formulae below are normalized by multiplying through by
467: local volume element (i.e. hx*hy) to obtain coefficients O(1) in two dimensions.
469:
470: */
471: dhx = (PetscReal)(info->mx-1); dhy = (PetscReal)(info->my-1);
472: hx = 1.0/dhx; hy = 1.0/dhy;
473: hxdhy = hx*dhy; /* hx/hy */ hydhx = hy*dhx; /* hy/hx */
475: xints = info->xs; xinte = info->xs+info->xm; yints = info->ys; yinte = info->ys+info->ym;
477: /* Test whether we are on the bottom edge of the global array */
478: if (yints == 0) {
479: j = 0;
480: yints = yints + 1;
481: /* bottom edge:
482: U = V = 0; Omega = -Grad_y(U)+Grad_x(V); dT/dn = 0; */
483: for (i=info->xs; i<info->xs+info->xm; i++) {
484: f[j][i].u = x[j][i].u;
485: f[j][i].v = x[j][i].v;
486: f[j][i].omega = x[j][i].omega + (x[j+1][i].u - x[j][i].u)*dhy;
487: f[j][i].temp = x[j][i].temp-x[j+1][i].temp;
488: }
489: }
491: /* Test whether we are on the top edge of the global array */
492: if (yinte == info->my) {
493: j = info->my - 1;
494: yinte = yinte - 1;
495: /* top edge:
496: U = lid; V = 0; Omega = -Grad_y(U)+Grad_x(V); dT/dn = 0; */
497: for (i=info->xs; i<info->xs+info->xm; i++) {
498: f[j][i].u = x[j][i].u - lid;
499: f[j][i].v = x[j][i].v;
500: f[j][i].omega = x[j][i].omega + (x[j][i].u - x[j-1][i].u)*dhy;
501: f[j][i].temp = x[j][i].temp-x[j-1][i].temp;
502: }
503: }
505: /* Test whether we are on the left edge of the global array */
506: if (xints == 0) {
507: i = 0;
508: xints = xints + 1;
509: /* left edge:
510: U = V = 0; Omega = -Grad_y(U)+Grad_x(V); T = 0; */
511: for (j=info->ys; j<info->ys+info->ym; j++) {
512: f[j][i].u = x[j][i].u;
513: f[j][i].v = x[j][i].v;
514: f[j][i].omega = x[j][i].omega - (x[j][i+1].v - x[j][i].v)*dhx;
515: f[j][i].temp = x[j][i].temp;
516: }
517: }
519: /* Test whether we are on the right edge of the global array */
520: if (xinte == info->mx) {
521: i = info->mx - 1;
522: xinte = xinte - 1;
523: /* right edge:
524: U = V = 0; Omega = -Grad_y(U)+Grad_x(V); T = grashof; */
525: for (j=info->ys; j<info->ys+info->ym; j++) {
526: f[j][i].u = x[j][i].u;
527: f[j][i].v = x[j][i].v;
528: f[j][i].omega = x[j][i].omega - (x[j][i].v - x[j][i-1].v)*dhx;
529: f[j][i].temp = x[j][i].temp - (PetscReal)(grashof>0);
530: }
531: }
533: /* Compute over the interior points */
534: for (j=yints; j<yinte; j++) {
535: for (i=xints; i<xinte; i++) {
536: /* convective coefficients for upwinding */
537: vx = x[j][i].u; avx = PetscAbsScalar(vx);
538: vxp = .5*(vx+avx); vxm = .5*(vx-avx);
539: vy = x[j][i].v; avy = PetscAbsScalar(vy);
540: vyp = .5*(vy+avy); vym = .5*(vy-avy);
542: /* U velocity */
543: u = x[j][i].u;
544: uxx = (2.0*u - x[j][i-1].u - x[j][i+1].u)*hydhx;
545: uyy = (2.0*u - x[j-1][i].u - x[j+1][i].u)*hxdhy;
546: f[j][i].u = uxx + uyy - .5*(x[j+1][i].omega-x[j-1][i].omega)*hx;
548: /* V velocity */
549: u = x[j][i].v;
550: uxx = (2.0*u - x[j][i-1].v - x[j][i+1].v)*hydhx;
551: uyy = (2.0*u - x[j-1][i].v - x[j+1][i].v)*hxdhy;
552: f[j][i].v = uxx + uyy + .5*(x[j][i+1].omega-x[j][i-1].omega)*hy;
554: /* Omega */
555: u = x[j][i].omega;
556: uxx = (2.0*u - x[j][i-1].omega - x[j][i+1].omega)*hydhx;
557: uyy = (2.0*u - x[j-1][i].omega - x[j+1][i].omega)*hxdhy;
558: f[j][i].omega = uxx + uyy
559: + (vxp*(u - x[j][i-1].omega) + vxm*(x[j][i+1].omega - u)) * hy
560: + (vyp*(u - x[j-1][i].omega) + vym*(x[j+1][i].omega - u)) * hx
561: - .5 * grashof * (x[j][i+1].temp - x[j][i-1].temp) * hy;
563: /* Temperature */
564: u = x[j][i].temp;
565: uxx = (2.0*u - x[j][i-1].temp - x[j][i+1].temp)*hydhx;
566: uyy = (2.0*u - x[j-1][i].temp - x[j+1][i].temp)*hxdhy;
567: f[j][i].temp = uxx + uyy
568: + prandtl * ((vxp*(u - x[j][i-1].temp) + vxm*(x[j][i+1].temp - u)) * hy
569: + (vyp*(u - x[j-1][i].temp) + vym*(x[j+1][i].temp - u)) * hx);
570: }
571: }
573: /* Flop count (multiply-adds are counted as 2 operations) */
574: PetscLogFlops(84.0*info->ym*info->xm);
575: return(0);
576: }
580: /*
581: Form function for Physics 1:
582: same as FormFunctionLocal() except without f.temp and x.temp.
583: the input x.temp comes from the solu_true
584: */
585: PetscErrorCode FormFunctionLocal1(DMDALocalInfo *info,Field1 **x,Field1 **f,void *ptr)
586: {
587: AppCtx *user = (AppCtx*)ptr;
588: PetscInt xints,xinte,yints,yinte,i,j;
589: PetscReal hx,hy,dhx,dhy,hxdhy,hydhx;
590: PetscReal grashof,lid; /* ,prandtl */
591: PetscScalar u,uxx,uyy,vx,vy,avx,avy,vxp,vxm,vyp,vym;
592: Field **solu=user->x;
593: Field2 **solu2=user->x2;
596: grashof = user->grashof;
597: /* prandtl = user->prandtl; */
598: lid = user->lidvelocity;
600: dhx = (PetscReal)(info->mx-1); dhy = (PetscReal)(info->my-1);
601: hx = 1.0/dhx; hy = 1.0/dhy;
602: hxdhy = hx*dhy; hydhx = hy*dhx;
604: xints = info->xs; xinte = info->xs+info->xm; yints = info->ys; yinte = info->ys+info->ym;
606: /* Test whether we are on the bottom edge of the global array */
607: if (yints == 0) {
608: j = 0;
609: yints = yints + 1;
610: /* bottom edge */
611: for (i=info->xs; i<info->xs+info->xm; i++) {
612: f[j][i].u = x[j][i].u;
613: f[j][i].v = x[j][i].v;
614: f[j][i].omega = x[j][i].omega + (x[j+1][i].u - x[j][i].u)*dhy;
615: }
616: }
618: /* Test whether we are on the top edge of the global array */
619: if (yinte == info->my) {
620: j = info->my - 1;
621: yinte = yinte - 1;
622: /* top edge */
623: for (i=info->xs; i<info->xs+info->xm; i++) {
624: f[j][i].u = x[j][i].u - lid;
625: f[j][i].v = x[j][i].v;
626: f[j][i].omega = x[j][i].omega + (x[j][i].u - x[j-1][i].u)*dhy;
627: }
628: }
630: /* Test whether we are on the left edge of the global array */
631: if (xints == 0) {
632: i = 0;
633: xints = xints + 1;
634: /* left edge */
635: for (j=info->ys; j<info->ys+info->ym; j++) {
636: f[j][i].u = x[j][i].u;
637: f[j][i].v = x[j][i].v;
638: f[j][i].omega = x[j][i].omega - (x[j][i+1].v - x[j][i].v)*dhx;
639: }
640: }
642: /* Test whether we are on the right edge of the global array */
643: if (xinte == info->mx) {
644: i = info->mx - 1;
645: xinte = xinte - 1;
646: /* right edge */
647: for (j=info->ys; j<info->ys+info->ym; j++) {
648: f[j][i].u = x[j][i].u;
649: f[j][i].v = x[j][i].v;
650: f[j][i].omega = x[j][i].omega - (x[j][i].v - x[j][i-1].v)*dhx;
651: }
652: }
654: /* Compute over the interior points */
655: for (j=yints; j<yinte; j++) {
656: for (i=xints; i<xinte; i++) {
657: /* convective coefficients for upwinding */
658: vx = x[j][i].u; avx = PetscAbsScalar(vx);
659: vxp = .5*(vx+avx); vxm = .5*(vx-avx);
660: vy = x[j][i].v; avy = PetscAbsScalar(vy);
661: vyp = .5*(vy+avy); vym = .5*(vy-avy);
663: /* U velocity */
664: u = x[j][i].u;
665: uxx = (2.0*u - x[j][i-1].u - x[j][i+1].u)*hydhx;
666: uyy = (2.0*u - x[j-1][i].u - x[j+1][i].u)*hxdhy;
667: f[j][i].u = uxx + uyy - .5*(x[j+1][i].omega-x[j-1][i].omega)*hx;
669: /* V velocity */
670: u = x[j][i].v;
671: uxx = (2.0*u - x[j][i-1].v - x[j][i+1].v)*hydhx;
672: uyy = (2.0*u - x[j-1][i].v - x[j+1][i].v)*hxdhy;
673: f[j][i].v = uxx + uyy + .5*(x[j][i+1].omega-x[j][i-1].omega)*hy;
675: /* Omega */
676: u = x[j][i].omega;
677: uxx = (2.0*u - x[j][i-1].omega - x[j][i+1].omega)*hydhx;
678: uyy = (2.0*u - x[j-1][i].omega - x[j+1][i].omega)*hxdhy;
679: f[j][i].omega = uxx + uyy
680: + (vxp*(u - x[j][i-1].omega) + vxm*(x[j][i+1].omega - u)) * hy
681: + (vyp*(u - x[j-1][i].omega) + vym*(x[j+1][i].omega - u)) * hx;
682: if (user->COMPOSITE_MODEL){
683: f[j][i].omega += - .5 * grashof * (solu2[j][i+1].temp - solu2[j][i-1].temp) * hy;
684: } else {
685: f[j][i].omega += - .5 * grashof * (solu[j][i+1].temp - solu[j][i-1].temp) * hy;
686: }
687: }
688: }
689: return(0);
690: }
694: /*
695: Form function for Physics 2:
696: same as FormFunctionLocal() but only has f.temp and x.temp.
697: the input x.u and x.v come from solu_true
698: */
699: PetscErrorCode FormFunctionLocal2(DMDALocalInfo *info,Field2 **x,Field2 **f,void *ptr)
700: {
701: AppCtx *user = (AppCtx*)ptr;
702: PetscInt xints,xinte,yints,yinte,i,j;
703: PetscReal hx,hy,dhx,dhy,hxdhy,hydhx;
704: PetscReal grashof,prandtl; /* ,lid */
705: PetscScalar u,uxx,uyy,vx,vy,avx,avy,vxp,vxm,vyp,vym;
706: Field **solu=user->x;
707: Field1 **solu1=user->x1;
710: grashof = user->grashof;
711: prandtl = user->prandtl;
712: /* lid = user->lidvelocity; */
714: dhx = (PetscReal)(info->mx-1); dhy = (PetscReal)(info->my-1);
715: hx = 1.0/dhx; hy = 1.0/dhy;
716: hxdhy = hx*dhy; hydhx = hy*dhx;
718: xints = info->xs; xinte = info->xs+info->xm; yints = info->ys; yinte = info->ys+info->ym;
720: /* Test whether we are on the bottom edge of the global array */
721: if (yints == 0) {
722: j = 0;
723: yints = yints + 1;
724: /* bottom edge */
725: for (i=info->xs; i<info->xs+info->xm; i++) {
726: f[j][i].temp = x[j][i].temp-x[j+1][i].temp;
727: }
728: }
730: /* Test whether we are on the top edge of the global array */
731: if (yinte == info->my) {
732: j = info->my - 1;
733: yinte = yinte - 1;
734: /* top edge */
735: for (i=info->xs; i<info->xs+info->xm; i++) {
736: f[j][i].temp = x[j][i].temp-x[j-1][i].temp;
737: }
738: }
740: /* Test whether we are on the left edge of the global array */
741: if (xints == 0) {
742: i = 0;
743: xints = xints + 1;
744: /* left edge */
745: for (j=info->ys; j<info->ys+info->ym; j++) {
746: f[j][i].temp = x[j][i].temp;
747: }
748: }
750: /* Test whether we are on the right edge of the global array */
751: if (xinte == info->mx) {
752: i = info->mx - 1;
753: xinte = xinte - 1;
754: /* right edge */
755: for (j=info->ys; j<info->ys+info->ym; j++) {
756: f[j][i].temp = x[j][i].temp - (PetscReal)(grashof>0);
757: }
758: }
760: /* Compute over the interior points */
761: for (j=yints; j<yinte; j++) {
762: for (i=xints; i<xinte; i++) {
763: /* convective coefficients for upwinding */
764: if (user->COMPOSITE_MODEL){
765: vx = solu1[j][i].u; vy = solu1[j][i].v;
766: } else {
767: vx = solu[j][i].u; vy = solu[j][i].v;
768: }
769: avx = PetscAbsScalar(vx);
770: vxp = .5*(vx+avx); vxm = .5*(vx-avx);
771: avy = PetscAbsScalar(vy);
772: vyp = .5*(vy+avy);
773: vym = .5*(vy-avy);
775: /* Temperature */
776: u = x[j][i].temp;
777: uxx = (2.0*u - x[j][i-1].temp - x[j][i+1].temp)*hydhx;
778: uyy = (2.0*u - x[j-1][i].temp - x[j+1][i].temp)*hxdhy;
779: f[j][i].temp = uxx + uyy
780: + prandtl * ((vxp*(u - x[j][i-1].temp) + vxm*(x[j][i+1].temp - u)) * hy
781: + (vyp*(u - x[j-1][i].temp) + vym*(x[j+1][i].temp - u)) * hx);
782: }
783: }
784: return(0);
785: }
789: /*
790: FormFunctionComp - Unwraps the input vector and passes its local ghosted pieces into the user function
791: */
792: PetscErrorCode FormFunctionComp(SNES snes,Vec X,Vec F,void *ctx)
793: {
795: DMMG dmmg = (DMMG)ctx;
796: AppCtx *user = (AppCtx*)dmmg->user;
797: DM dm = dmmg->dm;
798: DMDALocalInfo info1,info2;
799: DM da1,da2;
800: Field1 **x1,**f1;
801: Field2 **x2,**f2;
802: Vec X1,X2,F1,F2;
805: DMCompositeGetEntries(dm,&da1,&da2);
806: DMDAGetLocalInfo(da1,&info1);
807: DMDAGetLocalInfo(da2,&info2);
809: /* Get local vectors to hold ghosted parts of X */
810: DMCompositeGetLocalVectors(dm,&X1,&X2);
811: DMCompositeScatter(dm,X,X1,X2);
813: /* Access the arrays inside the subvectors of X */
814: DMDAVecGetArray(da1,X1,(void**)&x1);
815: DMDAVecGetArray(da2,X2,(void**)&x2);
817: /* Access the subvectors in F.
818: These are not ghosted and directly access the memory locations in F */
819: DMCompositeGetAccess(dm,F,&F1,&F2);
821: /* Access the arrays inside the subvectors of F */
822: DMDAVecGetArray(da1,F1,(void**)&f1);
823: DMDAVecGetArray(da2,F2,(void**)&f2);
825: /* Evaluate local user provided function */
826: user->x2 = x2; /* used by FormFunctionLocal1() */
827: FormFunctionLocal1(&info1,x1,f1,(void**)user);
828: user->x1 = x1; /* used by FormFunctionLocal2() */
829: FormFunctionLocal2(&info2,x2,f2,(void**)user);
831: DMDAVecRestoreArray(da1,F1,(void**)&f1);
832: DMDAVecRestoreArray(da2,F2,(void**)&f2);
833: DMCompositeRestoreAccess(dm,F,&F1,&F2);
834: DMDAVecRestoreArray(da1,X1,(void**)&x1);
835: DMDAVecRestoreArray(da2,X2,(void**)&x2);
836: DMCompositeRestoreLocalVectors(dm,&X1,&X2);
837: return(0);
838: }
842: /*
843: Visualize solutions
844: */
845: PetscErrorCode MySolutionView(MPI_Comm comm,PetscInt phy_num,void *ctx)
846: {
848: AppCtx *user = (AppCtx*)ctx;
849: DMMG *dmmg = user->dmmg;
850: DM da=DMMGGetDM(dmmg);
851: #if !defined(PETSC_USE_COMPLEX)
852: Field **x = user->x;
853: #endif
854: PetscInt i,j,xs,ys,xm,ym;
855: PetscMPIInt size;
858: #if defined(PETSC_USE_COMPLEX)
859: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"This routine does not support for scalar type complex yet\n");
860: #endif
861: MPI_Comm_size(PETSC_COMM_WORLD,&size);
862: DMDAGetCorners(da,&xs,&ys,PETSC_NULL,&xm,&ym,PETSC_NULL);
863:
864: switch (phy_num){
865: case 0:
866: if (size == 1){
867: PetscPrintf(PETSC_COMM_SELF,"Original Physics %d U,V,Omega,Temp: \n",phy_num);
868: PetscPrintf(PETSC_COMM_SELF,"-----------------------------------\n");
869: for (j=ys; j<ys+ym; j++) {
870: for (i=xs; i<xs+xm; i++) {
871: #if !defined(PETSC_USE_COMPLEX)
872: PetscPrintf(PETSC_COMM_SELF,"x[%d,%d] = %g, %g, %g, %g\n",j,i,x[j][i].u,x[j][i].v,x[j][i].omega,x[j][i].temp);
873: #endif
874: }
875: }
876: }
877: break;
878: case 1:
879: if (size == 1){
880: DMMG *dmmg1=user->dmmg1;
881: Vec solu_true = DMMGGetx(dmmg1);
882: DM da=DMMGGetDM(dmmg1);
883: Field1 **x1;
884: PetscPrintf(PETSC_COMM_SELF,"SubPhysics %d: U,V,Omega: \n",phy_num);
885: PetscPrintf(PETSC_COMM_SELF,"------------------------\n");
886: DMDAVecGetArray(da,solu_true,&x1);
887: for (j=ys; j<ys+ym; j++) {
888: for (i=xs; i<xs+xm; i++) {
889: #if !defined(PETSC_USE_COMPLEX)
890: PetscPrintf(PETSC_COMM_SELF,"x[%d,%d] = %g, %g, %g\n",j,i,x1[j][i].u,x1[j][i].v,x1[j][i].omega);
891: #endif
892: }
893: }
894: DMDAVecRestoreArray(da,solu_true,&x1);
895: }
896: break;
897: case 2:
898: if (size == 1){
899: DMMG *dmmg2=user->dmmg2;
900: Vec solu_true = DMMGGetx(dmmg2);
901: DM da=DMMGGetDM(dmmg2);
902: Field2 **x2;
903: PetscPrintf(PETSC_COMM_SELF,"SubPhysics %d: Temperature: \n",phy_num);
904: PetscPrintf(PETSC_COMM_SELF,"--------------------------\n");
905: DMDAVecGetArray(da,solu_true,&x2);
906: for (j=ys; j<ys+ym; j++) {
907: for (i=xs; i<xs+xm; i++) {
908: #if !defined(PETSC_USE_COMPLEX)
909: PetscPrintf(PETSC_COMM_SELF,"x[%d,%d] = %g\n",j,i,x2[j][i].temp);
910: #endif
911: }
912: }
913: DMDAVecRestoreArray(da,solu_true,&x2);
914: }
915: break;
916: default:
917: if (size == 1){
918: DMMG *dmmg_comp=user->dmmg_comp;
919: DM da1,da2,da=DMMGGetDM(dmmg);
920: Vec X1,X2,solu_true = DMMGGetx(dmmg);
921: Field **x;
922: Field1 **x1;
923: Field2 **x2;
924: DM dm = (*dmmg_comp)->dm;
925: PetscReal err,err_tmp;
926: if (phy_num == 3){
927: PetscPrintf(PETSC_COMM_SELF,"Composite physics %d, U,V,Omega,Temp: \n",phy_num);
928: PetscPrintf(PETSC_COMM_SELF,"------------------------------------\n");
929: PetscPrintf(PETSC_COMM_SELF,"Composite physics, U,V,Omega,Temp: \n");
930: }
931: DMDAVecGetArray(da,solu_true,&x);
932: DMCompositeGetEntries(dm,&da1,&da2);
933: DMCompositeGetLocalVectors(dm,&X1,&X2);
934: DMDAVecGetArray(da1,X1,(void**)&x1);
935: DMDAVecGetArray(da2,X2,(void**)&x2);
937: err = 0.0;
938: for (j=ys; j<ys+ym; j++) {
939: for (i=xs; i<xs+xm; i++) {
940: err_tmp = PetscAbsScalar(x[j][i].u-x1[j][i].u) + PetscAbsScalar(x[j][i].v-x1[j][i].v) + PetscAbsScalar(x[j][i].omega-x1[j][i].omega);
941: err_tmp += PetscAbsScalar(x[j][i].temp-x2[j][i].temp);
942: if (err < err_tmp) err = err_tmp;
943: if (phy_num == 3){
944: #if !defined(PETSC_USE_COMPLEX)
945: PetscPrintf(PETSC_COMM_SELF,"x[%d,%d] = %g, %g, %g, %g\n",j,i,x1[j][i].u,x1[j][i].v,x1[j][i].omega,x2[j][i].temp);
946: #endif
947: }
948: }
949: }
950: PetscPrintf(PETSC_COMM_SELF,"|solu - solu_comp| = %g\n",err);
951: DMDAVecRestoreArray(da1,X1,(void**)&x1);
952: DMDAVecRestoreArray(da2,X2,(void**)&x2);
953: DMCompositeRestoreLocalVectors(dm,&X1,&X2);
954: DMDAVecRestoreArray(da,solu_true,&x);
955: }
956: }
957: return(0);
958: }