Actual source code: ex4tu.c
1: /* Program usage: mpiexec -n <procs> ex5 [-help] [all PETSc options] */
3: static char help[] = "Nonlinear PDE in 2d.\n\
4: We solve the Lane-Emden equation in a 2D rectangular\n\
5: domain, using distributed arrays (DMDAs) to partition the parallel grid.\n\n";
7: /*T
8: Concepts: SNES^parallel Lane-Emden example
9: Concepts: DMDA^using distributed arrays;
10: Processors: n
11: T*/
13: /* ------------------------------------------------------------------------
15: The Lane-Emden equation is given by the partial differential equation
16:
17: -alpha*Laplacian u - lambda*u^3 = 0, 0 < x,y < 1,
18:
19: with boundary conditions
20:
21: u = 0 for x = 0, x = 1, y = 0, y = 1.
22:
23: A bilinear finite element approximation is used to discretize the boundary
24: value problem to obtain a nonlinear system of equations.
26: ------------------------------------------------------------------------- */
28: /*
29: Include "petscdmda.h" so that we can use distributed arrays (DMDAs).
30: Include "petscsnes.h" so that we can use SNES solvers. Note that this
31: file automatically includes:
32: petscsys.h - base PETSc routines petscvec.h - vectors
33: petscmat.h - matrices
34: petscis.h - index sets petscksp.h - Krylov subspace methods
35: petscviewer.h - viewers petscpc.h - preconditioners
36: petscksp.h - linear solvers
37: */
38: #include <petscdmda.h>
39: #include <petscdmmg.h>
40: #include <petscsnes.h>
41: #include <../src/snes/impls/ls/lsimpl.h>
42: /*
43: User-defined application context - contains data needed by the
44: application-provided call-back routines, FormJacobianLocal() and
45: FormFunctionLocal().
46: */
47: typedef struct {
48: // DM da; /* distributed array data structure */
49: PetscReal alpha; /* parameter controlling linearity */
50: PetscReal lambda; /* parameter controlling nonlinearity */
51: PetscBool draw_contours; /* flag - 1 indicates drawing contours */
52: } AppCtx;
55: /*
56: User-defined routines
57: */
66: int main(int argc,char **argv)
67: {
68: DMMG *dmmg; /* multilevel grid structure */
69: SNES snes; /* nonlinear solver */
70: //Vec x,r; /* solution, residual vectors */
71: //Mat J; /* Jacobian matrix */
72: AppCtx user; /* user-defined work context */
73: PetscInt its; /* iterations for convergence */
74: SNESConvergedReason reason;
75: PetscErrorCode ierr;
76: PetscReal lambda_max = 6.81, lambda_min = 0.0;
77: MPI_Comm comm;
78: PetscInt mx,my;
79: DM da;
80: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
81: Initialize program
82: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
83: PetscInitialize(&argc,&argv,(char *)0,help);
84: comm = PETSC_COMM_WORLD;
86: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
87: Initialize problem parameters
88: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
89: user.alpha = 1.0;
90: user.lambda = 6.0;
91: PetscOptionsGetReal(PETSC_NULL,"-alpha",&user.alpha,PETSC_NULL);
92: PetscOptionsGetReal(PETSC_NULL,"-lambda",&user.lambda,PETSC_NULL);
93: if (user.lambda > lambda_max || user.lambda < lambda_min) {
94: SETERRQ3(PETSC_COMM_SELF,1,"Lambda %g is out of range [%g, %g]", user.lambda, lambda_min, lambda_max);
95: }
98: // in order only run once, I block it: PetscPreLoadBegin(PETSC_TRUE,"SetUp");
99: DMMGCreate(comm,2,&user,&dmmg);
102: /*
103: Create distributed array multigrid object (DMMG) to manage parallel grid and vectors
104: for principal unknowns (x) and governing residuals (f)
105: */
106: DMDACreate2d(PETSC_COMM_WORLD, DMDA_BOUNDARY_NONE, DMDA_BOUNDARY_NONE,DMDA_STENCIL_BOX,-4,-4,PETSC_DECIDE,PETSC_DECIDE,1,1,PETSC_NULL,PETSC_NULL,&da);
107:
108: DMMGSetDM(dmmg,(DM)da);
109: DMDestroy(&da);
111: DMDAGetInfo(DMMGGetDM(dmmg),0,&mx,&my,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,
112: PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE);
113: PetscPrintf(comm,"mx = %d, my= %d\n",
114: mx,my);
115:
116: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
117: Create user context, set problem data, create vector data structures.
118: Also, compute the initial guess.
119: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
121: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
122: Create nonlinear solver context
124:
125: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
126: DMMGSetSNESLocal(dmmg,FormFunctionLocal,0,ad_FormFunctionLocal,admf_FormFunctionLocal);
127: DMMGSetFromOptions(dmmg);
128: DMMGSetSNESLocali(dmmg,FormFunctionLocali,ad_FormFunctionLocali,admf_FormFunctionLocali);
129: DMMGSetSNESLocalib(dmmg,FormFunctionLocali4,ad_FormFunctionLocali4,admf_FormFunctionLocali4);
131:
133: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
134: Solve the nonlinear system
135: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
136: DMMGSetInitialGuess(dmmg,FormInitialGuess);
138: //I block it: PetscPreLoadStage("Solve");
139: DMMGSolve(dmmg);
141: snes = DMMGGetSNES(dmmg);
142: SNESGetIterationNumber(snes,&its);
143: PetscPrintf(comm,"Number of Newton iterations = %D\n", its);
144:
145: /*
146: Visualize solution
147: */
148: PetscOptionsHasName(PETSC_NULL,"-contours",&user.draw_contours);
149: if (user.draw_contours) {
150: VecView(DMMGGetx(dmmg),PETSC_VIEWER_DRAW_WORLD);
151: //VecView(DMMGGetx(dmmg),PETSC_VIEWER_STDOUT_WORLD);
152: }
154: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
155: Free work space. All PETSc objects should be destroyed when they
156: are no longer needed.
157: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
159: DMMGDestroy(dmmg);
160: //PetscPreLoadEnd();
162: PetscFinalize();
163: return 0;
165: }
169: /*
170: FormInitialGuess - Forms initial approximation.
172: Input Parameters:
173: user - user-defined application context
174: X - vector
176: Output Parameter:
177: X - vector
178: */
179: PetscErrorCode FormInitialGuess(DMMG dmmg,Vec X)
180: {
181: AppCtx *user = (AppCtx*)dmmg->user;
182: DM da = dmmg->dm;
183: PetscInt i,j,Mx,My,xs,ys,xm,ym;
185: PetscReal lambda,temp1,temp,hx,hy;
186: PetscScalar **x;
189: DMDAGetInfo(da,PETSC_IGNORE,&Mx,&My,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,
190: PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE);
191:
192: lambda = user->lambda;
193: hx = 1.0/(PetscReal)(Mx-1);
194: hy = 1.0/(PetscReal)(My-1);
195: if (lambda == 0.0) {
196: temp1 = 1.0;
197: } else {
198: temp1 = lambda/(lambda + 1.0);
199: }
201: /*
202: Get a pointer to vector data.
203: - For default PETSc vectors, VecGetArray() returns a pointer to
204: the data array. Otherwise, the routine is implementation dependent.
205: - You MUST call VecRestoreArray() when you no longer need access to
206: the array.
207: */
208: DMDAVecGetArray(da,X,&x);
210: /*
211: Get local grid boundaries (for 2-dimensional DMDA):
212: xs, ys - starting grid indices (no ghost points)
213: xm, ym - widths of local grid (no ghost points)
215: */
216: DMDAGetCorners(da,&xs,&ys,PETSC_NULL,&xm,&ym,PETSC_NULL);
218: /*
219: Compute initial guess over the locally owned part of the grid
220: */
221: for (j=ys; j<ys+ym; j++) {
222: temp = (PetscReal)(PetscMin(j,My-j-1))*hy;
223: for (i=xs; i<xs+xm; i++) {
225: if (i == 0 || j == 0 || i == Mx-1 || j == My-1) {
226: /* boundary conditions are all zero Dirichlet */
227: x[j][i] = 0.0;
228: } else {
229: x[j][i] = temp1*sqrt(PetscMin((PetscReal)(PetscMin(i,Mx-i-1))*hx,temp));
230: }
231: }
232: }
234: /*
235: Restore vector
236: */
237: DMDAVecRestoreArray(da,X,&x);
239: return(0);
240: }
245: /*
246: FormFunctionLocal - Evaluates nonlinear function, F(x).
248: Process adiC(36): FormFunctionLocal nonlinearResidual FormFunctionLocali FormFunctionLocali4
250: */
251: PetscErrorCode FormFunctionLocal(DMDALocalInfo *info,PetscScalar **x,PetscScalar **f,AppCtx *user)
253: {
254: PetscScalar uLocal[4];
255: PetscScalar rLocal[4];
256: PetscScalar lintx[4],linty[4],lintw[4],intnum=4, int1,int2,int3,int4,p1,p2,p3,p4,p;
257: PetscReal alpha,lambda,hx,hy,hxhy,sc;
258: PetscInt i,j,k,l;
261: /* Compute function over the locally owned part of the grid. For each
262: vertex (i,j), we consider the element below:
264: 3 2
265: i,j+1 --- i+1,j
266: | |
267: | |
268: i,j --- i+1,j
269: 0 1
271: and therefore we do not loop over the last vertex in each dimension.
272: */
274: alpha = user->alpha;
275: lambda = user->lambda;
276: hx = 1.0/(PetscReal)(info->mx-1);
277: hy = 1.0/(PetscReal)(info->my-1);
278: sc = hx*hy*lambda;
279: hxhy = hx*hy;
282: // set all function values to be zero, this maybe not good for parallel computing
283: for(j = info->ys; j < info->ys+info->ym-1; j++) {
284: for(i = info->xs; i < info->xs+info->xm-1; i++) {
285: f[j][i]=0.0;
286: }
287: }
289:
290: lintx[0] = 0.21132486540518713;//(1.0-sqrt(1.0/3.0))/2.0;
291: lintx[1] = 0.78867513459481287;//(1.0+sqrt(1.0/3.0))/2.0;
292: lintx[2] = 0.78867513459481287;//(1.0+sqrt(1.0/3.0))/2.0;
293: lintx[3] = 0.21132486540518713;//(1.0-sqrt(1.0/3.0))/2.0;
294: linty[0] = 0.21132486540518713;//(1.0-sqrt(1.0/3.0))/2.0;
295: linty[1] = 0.21132486540518713;//(1.0-sqrt(1.0/3.0))/2.0;
296: linty[2] = 0.78867513459481287;//(1.0+sqrt(1.0/3.0))/2.0;
297: linty[3] = 0.78867513459481287;//(1.0+sqrt(1.0/3.0))/2.0;
299: lintw[0] = 1.0/4.0;
300: lintw[1] = 1.0/4.0;
301: lintw[2] = 1.0/4.0;
302: lintw[3] = 1.0/4.0;
303:
304:
305: for(j = info->ys; j < info->ys+info->ym-1; j++) {
306: for(i = info->xs; i < info->xs+info->xm-1; i++) {
307: uLocal[0] = x[j][i];
308: uLocal[1] = x[j][i+1];
309: uLocal[2] = x[j+1][i+1];
310: uLocal[3] = x[j+1][i];
311:
312:
313: /* Laplace term */
314: rLocal[0] = 2.0/3.0*uLocal[0]-1.0/6.0*uLocal[1]-1.0/3.0*uLocal[2]-1.0/6.0*uLocal[3];
315: rLocal[0] *= hxhy*alpha;
317: rLocal[1] = -1.0/6.0*uLocal[0]+2.0/3.0*uLocal[1]-1.0/6.0*uLocal[2]-1.0/3.0*uLocal[3];
318: rLocal[1] *= hxhy*alpha;
320: rLocal[2] = -1.0/3.0*uLocal[0]-1.0/6.0*uLocal[1]+2.0/3.0*uLocal[2]-1.0/6.0*uLocal[3];
321: rLocal[2] *= hxhy*alpha;
323: rLocal[3] = -1.0/6.0*uLocal[0]-1.0/3.0*uLocal[1]-1.0/6.0*uLocal[2]+2.0/3.0*uLocal[3];
324: rLocal[3] *= hxhy*alpha;
325:
326: /* nonlinear term */
327: int1 = 0; int2 = 0; int3 = 0; int4 = 0;
328: for(k=0; k<intnum; k++){
329: p1 = (1.0-lintx[k])*(1.0-linty[k]);
330: p2 = lintx[k] *(1.0-linty[k]);
331: p3 = lintx[k] * linty[k];
332: p4 = (1.0-lintx[k])* linty[k];
333: p = lintw[k]*PetscExpScalar( uLocal[0]*p1+ uLocal[1]*p2 + uLocal[2]*p3+ uLocal[3]*p4 );
334: // p = lintw[k]*( uLocal[0]*p1+ uLocal[1]*p2 + uLocal[2]*p3+ uLocal[3]*p4 )*( uLocal[0]*p1+ uLocal[1]*p2 + uLocal[2]*p3+ uLocal[3]*p4 )*( uLocal[0]*p1+ uLocal[1]*p2 + uLocal[2]*p3+ uLocal[3]*p4 );
335: int1 = int1 + p*p1;
336: int2 = int2 + p*p2;
337: int3 = int3 + p*p3;
338: int4 = int4 + p*p4;
339: }
340:
341: f[j][i] += rLocal[0] + int1*hxhy*(-1.0*sc);
342: f[j][i+1] += rLocal[1] + int2*hxhy*(-1.0*sc);
343: f[j+1][i+1] += rLocal[2] + int3*hxhy*(-1.0*sc);
344: f[j+1][i] += rLocal[3] + int4*hxhy*(-1.0*sc);
345:
346: if (i == 0 || j == 0) {
347: f[j][i] = x[j][i];
348: }
349: if ((i == info->mx-2) || (j == 0)) {
350: f[j][i+1] = x[j][i+1];
351: }
352: if ((i == info->mx-2) || (j == info->my-2)) {
353: f[j+1][i+1] = x[j+1][i+1];
354: }
355: if ((i == 0) || (j == info->my-2)) {
356: f[j+1][i] = x[j+1][i];
357: }
358: }
359: }
361: PetscLogFlops(68.0*(info->ym-1)*(info->xm-1));
362: return(0);
363: }
368: PetscErrorCode FormFunctionLocali(DMDALocalInfo *info,MatStencil *st,PetscScalar **x,PetscScalar *f,AppCtx *user)
369: {
370: PetscScalar uLocal[4];
371: PetscScalar rLocal[4];
372: PetscReal alpha,lambda,hx,hy,hxhy,sc;
373: PetscScalar lintx[4],linty[4],lintw[4],intnum=4, int1,int2,int3,int4,p1,p2,p3,p4,p;
374: PetscInt i,j,k,l;
378: alpha = user->alpha;
379: lambda = user->lambda;
380: hx = 1.0/(PetscReal)(info->mx-1);
381: hy = 1.0/(PetscReal)(info->my-1);
382: sc = -1.0*hx*hy*lambda;
383: hxhy = hx*hy;
385: /* Compute function over the locally owned part of the grid. For each
386: vertex (i,j), we consider the element below:
388: 3 2
389: i,j+1 --- i+1,j
390: | |
391: | |
392: i,j --- i+1,j
393: 0 1
395: and therefore we do not loop over the last vertex in each dimension.
396: */
399:
400: lintx[0] = 0.21132486540518713;//(1.0-sqrt(1.0/3.0))/2.0;
401: lintx[1] = 0.78867513459481287;//(1.0+sqrt(1.0/3.0))/2.0;
402: lintx[2] = 0.78867513459481287;//(1.0+sqrt(1.0/3.0))/2.0;
403: lintx[3] = 0.21132486540518713;//(1.0-sqrt(1.0/3.0))/2.0;
404: linty[0] = 0.21132486540518713;//(1.0-sqrt(1.0/3.0))/2.0;
405: linty[1] = 0.21132486540518713;//(1.0-sqrt(1.0/3.0))/2.0;
406: linty[2] = 0.78867513459481287;//(1.0+sqrt(1.0/3.0))/2.0;
407: linty[3] = 0.78867513459481287;//(1.0+sqrt(1.0/3.0))/2.0;
409: lintw[0] = 1.0/4.0;
410: lintw[1] = 1.0/4.0;
411: lintw[2] = 1.0/4.0;
412: lintw[3] = 1.0/4.0;
413:
415: i = st->i; j = st->j;
417: /*
418:
419: i-1,j+1 --- i,j+1 --- i+1,j+1
420: | | |
421: | 2 (1) | 1(0) |
422: i-1,j --- i,j --- i+1,j
423: | | |
424: | 3(2) | 4(3) |
425: i-1,j-1 --- i,j-1--- i+1,j-1
426:
427: */
428:
430: // boundary treatment
432: if (i == 0 || j == 0 || i == info->mx-1 || j== info->my-1) {
433: *f = x[j][i];
434: return(0);
435: }
436:
438: /* element 1 */
439: uLocal[0] = x[j][i];
440: uLocal[1] = x[j][i+1];
441: uLocal[2] = x[j+1][i+1];
442: uLocal[3] = x[j+1][i];
443:
444: rLocal[0] = 2.0/3.0*uLocal[0]-1.0/6.0*uLocal[1]-1.0/3.0*uLocal[2]-1.0/6.0*uLocal[3];
445: rLocal[0] *= hxhy*alpha;
446:
447:
448: int1 = 0;
449: for(k=0; k<intnum; k++){
450: p1 = (1.0-lintx[k])*(1.0-linty[k]);
451: p2 = lintx[k] *(1.0-linty[k]);
452: p3 = lintx[k] * linty[k];
453: p4 = (1.0-lintx[k])* linty[k];
454: p = lintw[k]*PetscExpScalar( uLocal[0]*p1+ uLocal[1]*p2 + uLocal[2]*p3+ uLocal[3]*p4 );
455: int1 = int1 + p*p1;
456: }
457:
458: *f = rLocal[0] + int1*hxhy*(sc);
460: /* element 2 */
461: uLocal[0] = x[j][i-1];
462: uLocal[1] = x[j][i];
463: uLocal[2] = x[j+1][i];
464: uLocal[3] = x[j+1][i-1];
465:
466: rLocal[1] = -1.0/6.0*uLocal[0]+2.0/3.0*uLocal[1]-1.0/6.0*uLocal[2]-1.0/3.0*uLocal[3];
467: rLocal[1] *= hxhy*alpha;
468: int2 = 0;
469: for(k=0; k<intnum; k++){
470: p1 = (1.0-lintx[k])*(1.0-linty[k]);
471: p2 = lintx[k] *(1.0-linty[k]);
472: p3 = lintx[k] * linty[k];
473: p4 = (1.0-lintx[k])* linty[k];
474: p = lintw[k]*PetscExpScalar( uLocal[0]*p1+ uLocal[1]*p2 + uLocal[2]*p3+ uLocal[3]*p4 );
475: int2 = int2 + p*p2;
476: }
477:
478: *f += rLocal[1] + int2*hxhy*(sc);
479:
480: /* element 3 */
481: uLocal[0] = x[j-1][i-1];
482: uLocal[1] = x[j-1][i];
483: uLocal[2] = x[j][i];
484: uLocal[3] = x[j][i-1];
485:
486: rLocal[2] = -1.0/3.0*uLocal[0]-1.0/6.0*uLocal[1]+2.0/3.0*uLocal[2]-1.0/6.0*uLocal[3];
487: rLocal[2] *= hxhy*alpha;
488: int3 = 0;
489: for(k=0; k<intnum; k++){
490: p1 = (1.0-lintx[k])*(1.0-linty[k]);
491: p2 = lintx[k] *(1.0-linty[k]);
492: p3 = lintx[k] * linty[k];
493: p4 = (1.0-lintx[k])* linty[k];
494: p = lintw[k]*PetscExpScalar( uLocal[0]*p1+ uLocal[1]*p2 + uLocal[2]*p3+ uLocal[3]*p4 );
495: int3 = int3 + p*p3;
496: }
497:
498: *f += rLocal[2] + int3*hxhy*(sc);
500: /* element 4 */
501: uLocal[0] = x[j-1][i];
502: uLocal[1] = x[j-1][i+1];
503: uLocal[2] = x[j][i+1];
504: uLocal[3] = x[j][i];
506: rLocal[3] = -1.0/6.0*uLocal[0]-1.0/3.0*uLocal[1]-1.0/6.0*uLocal[2]+2.0/3.0*uLocal[3];
507: rLocal[3] *= hxhy*alpha;
508: int4 = 0;
509: for(k=0; k<intnum; k++){
510: p1 = (1.0-lintx[k])*(1.0-linty[k]);
511: p2 = lintx[k] *(1.0-linty[k]);
512: p3 = lintx[k] * linty[k];
513: p4 = (1.0-lintx[k])* linty[k];
514: p = lintw[k]*PetscExpScalar( uLocal[0]*p1+ uLocal[1]*p2 + uLocal[2]*p3+ uLocal[3]*p4 );
515: int4 = int4 + p*p4;
516: }
517:
518: *f += rLocal[3] + int4*hxhy*(sc);
522: return(0);
523: }
529: //****************************************************************************
535: PetscErrorCode FormFunctionLocali4(DMDALocalInfo *info,MatStencil *st,PetscScalar **x,PetscScalar *f,AppCtx *user)
536: {
537: PetscScalar uLocal[4],fLocal[5][5];
538: PetscScalar rLocal[4];
539: PetscScalar lintx[4],linty[4],lintw[4],intnum=4, int1,int2,int3,int4,p1,p2,p3,p4,p;
540: PetscReal alpha,lambda,hx,hy,hxhy,sc;
541: PetscInt i,j,k,kk,ll,istar,iend,jstar,jend,imax,jmax,id,jd;
546: alpha = user->alpha;
547: lambda = user->lambda;
548: hx = 1.0/(PetscReal)(info->mx-1);
549: hy = 1.0/(PetscReal)(info->my-1);
550: sc = -1.0*hx*hy*lambda;
551: hxhy = hx*hy;
553: /* Compute function over the locally owned part of the grid. For each
554: vertex (i,j), we consider the element below:
556: 3 2
557: i,j+1 --- i+1,j
558: | |
559: | |
560: i,j --- i+1,j
561: 0 1
563: and therefore we do not loop over the last vertex in each dimension.
564: */
567:
569: i = st->i; j = st->j;
571: /*
572:
573: i-1,j+1 --- i,j+1 --- i+1,j+1
574: | | |
575: | 2 (1) | 1(0) |
576: i-1,j --- i,j --- i+1,j
577: | | |
578: | 3(2) | 4(3) |
579: i-1,j-1 --- i,j-1--- i+1,j-1
580:
581: */
582:
584: // boundary treatment
586: if (i == 0 || j == 0 || i == info->mx-1 || j== info->my-1) {
587: f[0] = x[j][i];
588: return(0);
589: }
590:
591: //(2) second cases: next to boundary
594: if(i==1 || j==1 ||i==info->mx-2 ||j==info->my-2){
595:
597: lintx[0] = 0.21132486540518713;//(1.0-sqrt(1.0/3.0))/2.0;
598: lintx[1] = 0.78867513459481287;//(1.0+sqrt(1.0/3.0))/2.0;
599: lintx[2] = 0.78867513459481287;//(1.0+sqrt(1.0/3.0))/2.0;
600: lintx[3] = 0.21132486540518713;//(1.0-sqrt(1.0/3.0))/2.0;
601: linty[0] = 0.21132486540518713;//(1.0-sqrt(1.0/3.0))/2.0;
602: linty[1] = 0.21132486540518713;//(1.0-sqrt(1.0/3.0))/2.0;
603: linty[2] = 0.78867513459481287;//(1.0+sqrt(1.0/3.0))/2.0;
604: linty[3] = 0.78867513459481287;//(1.0+sqrt(1.0/3.0))/2.0;
606: lintw[0] = 1.0/4.0;
607: lintw[1] = 1.0/4.0;
608: lintw[2] = 1.0/4.0;
609: lintw[3] = 1.0/4.0;
610:
612: /* element 1 */
613: uLocal[0] = x[j][i];
614: uLocal[1] = x[j][i+1];
615: uLocal[2] = x[j+1][i+1];
616: uLocal[3] = x[j+1][i];
617:
618: rLocal[0] = 2.0/3.0*uLocal[0]-1.0/6.0*uLocal[1]-1.0/3.0*uLocal[2]-1.0/6.0*uLocal[3];
619: rLocal[0] *= hxhy*alpha;
620:
621:
622: int1 = 0;
623: for(k=0; k<intnum; k++){
624: p1 = (1.0-lintx[k])*(1.0-linty[k]);
625: p2 = lintx[k] *(1.0-linty[k]);
626: p3 = lintx[k] * linty[k];
627: p4 = (1.0-lintx[k])* linty[k];
628: p = lintw[k]*PetscExpScalar( uLocal[0]*p1+ uLocal[1]*p2 + uLocal[2]*p3+ uLocal[3]*p4 );
629: int1 = int1 + p*p1;
630: }
631:
632: f[0] = rLocal[0] + int1*hxhy*(sc);
634: /* element 2 */
635: uLocal[0] = x[j][i-1];
636: uLocal[1] = x[j][i];
637: uLocal[2] = x[j+1][i];
638: uLocal[3] = x[j+1][i-1];
639:
640: rLocal[1] = -1.0/6.0*uLocal[0]+2.0/3.0*uLocal[1]-1.0/6.0*uLocal[2]-1.0/3.0*uLocal[3];
641: rLocal[1] *= hxhy*alpha;
642: int2 = 0;
643: for(k=0; k<intnum; k++){
644: p1 = (1.0-lintx[k])*(1.0-linty[k]);
645: p2 = lintx[k] *(1.0-linty[k]);
646: p3 = lintx[k] * linty[k];
647: p4 = (1.0-lintx[k])* linty[k];
648: p = lintw[k]*PetscExpScalar( uLocal[0]*p1+ uLocal[1]*p2 + uLocal[2]*p3+ uLocal[3]*p4 );
649: int2 = int2 + p*p2;
650: }
651:
652: f[0] += rLocal[1] + int2*hxhy*(sc);
653:
654: /* element 3 */
655: uLocal[0] = x[j-1][i-1];
656: uLocal[1] = x[j-1][i];
657: uLocal[2] = x[j][i];
658: uLocal[3] = x[j][i-1];
659:
660: rLocal[2] = -1.0/3.0*uLocal[0]-1.0/6.0*uLocal[1]+2.0/3.0*uLocal[2]-1.0/6.0*uLocal[3];
661: rLocal[2] *= hxhy*alpha;
662: int3 = 0;
663: for(k=0; k<intnum; k++){
664: p1 = (1.0-lintx[k])*(1.0-linty[k]);
665: p2 = lintx[k] *(1.0-linty[k]);
666: p3 = lintx[k] * linty[k];
667: p4 = (1.0-lintx[k])* linty[k];
668: p = lintw[k]*PetscExpScalar( uLocal[0]*p1+ uLocal[1]*p2 + uLocal[2]*p3+ uLocal[3]*p4 );
669: int3 = int3 + p*p3;
670: }
671:
672: f[0] += rLocal[2] + int3*hxhy*(sc);
674: /* element 4 */
675: uLocal[0] = x[j-1][i];
676: uLocal[1] = x[j-1][i+1];
677: uLocal[2] = x[j][i+1];
678: uLocal[3] = x[j][i];
680: rLocal[3] = -1.0/6.0*uLocal[0]-1.0/3.0*uLocal[1]-1.0/6.0*uLocal[2]+2.0/3.0*uLocal[3];
681: rLocal[3] *= hxhy*alpha;
682: int4 = 0;
683: for(k=0; k<intnum; k++){
684: p1 = (1.0-lintx[k])*(1.0-linty[k]);
685: p2 = lintx[k] *(1.0-linty[k]);
686: p3 = lintx[k] * linty[k];
687: p4 = (1.0-lintx[k])* linty[k];
688: p = lintw[k]*PetscExpScalar( uLocal[0]*p1+ uLocal[1]*p2 + uLocal[2]*p3+ uLocal[3]*p4 );
689: int4 = int4 + p*p4;
690: }
691:
692: f[0] += rLocal[3] + int4*hxhy*(sc);
696: return(0);
697: }
699: /*
700: i-2,j+2 --- i-1,j+2 --- i,j+2 --- i+1,j+2 --- i+2,j+2
701: | | | | |
702: | e13 | e14 | e15 | e16 |
703: i-2,j+1 --- i-1,j+1 --- i,j+1 --- i+1,j+1 --- i+2;j+1
704: | |(6) |(7) |(8) |
705: | e9 | e10 | e11 | e12 |
706: i-2,j --- i-1,j --- i,j --- i+1,j --- i+2,j
707: | |(3) |(4) |(5) |
708: | e5 | e6 | e7 | e8 |
709: i-2,j-1 --- i-1,j-1 --- i,j-1--- i+1,j-1 --- i+2,j-1
710: | |(0) |(1) |(2) |
711: | e1 | e2 | e3 | e4 |
712: i-2,j-2 --- i-1,j-2 --- i,j-2--- i+1,j-2 --- i+2,j-2
713:
714: */
715:
716: jstar = j - 2;
717: jend = j + 2;
718: istar = i - 2;
719: iend = i + 2;
720:
721:
722: lintx[0] = 0.21132486540518713;//(1.0-sqrt(1.0/3.0))/2.0;
723: lintx[1] = 0.78867513459481287;//(1.0+sqrt(1.0/3.0))/2.0;
724: lintx[2] = 0.78867513459481287;//(1.0+sqrt(1.0/3.0))/2.0;
725: lintx[3] = 0.21132486540518713;//(1.0-sqrt(1.0/3.0))/2.0;
726: linty[0] = 0.21132486540518713;//(1.0-sqrt(1.0/3.0))/2.0;
727: linty[1] = 0.21132486540518713;//(1.0-sqrt(1.0/3.0))/2.0;
728: linty[2] = 0.78867513459481287;//(1.0+sqrt(1.0/3.0))/2.0;
729: linty[3] = 0.78867513459481287;//(1.0+sqrt(1.0/3.0))/2.0;
731: lintw[0] = 1.0/4.0;
732: lintw[1] = 1.0/4.0;
733: lintw[2] = 1.0/4.0;
734: lintw[3] = 1.0/4.0;
735:
737: id=0;
738: for(kk=1;kk<4;kk++) {
739: for(ll=1; ll<4;ll++) {
740: i=istar+ll;
741: j=jstar+kk;
743: /* element 1 */
744: uLocal[0] = x[j][i];
745: uLocal[1] = x[j][i+1];
746: uLocal[2] = x[j+1][i+1];
747: uLocal[3] = x[j+1][i];
748:
749: rLocal[0] = 2.0/3.0*uLocal[0]-1.0/6.0*uLocal[1]-1.0/3.0*uLocal[2]-1.0/6.0*uLocal[3];
750: rLocal[0] *= hxhy*alpha;
751:
752:
753: int1 = 0;
754: for(k=0; k<intnum; k++){
755: p1 = (1.0-lintx[k])*(1.0-linty[k]);
756: p2 = lintx[k] *(1.0-linty[k]);
757: p3 = lintx[k] * linty[k];
758: p4 = (1.0-lintx[k])* linty[k];
759: p = lintw[k]*PetscExpScalar( uLocal[0]*p1+ uLocal[1]*p2 + uLocal[2]*p3+ uLocal[3]*p4 );
760: int1 = int1 + p*p1;
761: }
762:
763: f[id] = rLocal[0] + int1*hxhy*(sc);
765: /* element 2 */
766: uLocal[0] = x[j][i-1];
767: uLocal[1] = x[j][i];
768: uLocal[2] = x[j+1][i];
769: uLocal[3] = x[j+1][i-1];
770:
771: rLocal[1] = -1.0/6.0*uLocal[0]+2.0/3.0*uLocal[1]-1.0/6.0*uLocal[2]-1.0/3.0*uLocal[3];
772: rLocal[1] *= hxhy*alpha;
773: int2 = 0;
774: for(k=0; k<intnum; k++){
775: p1 = (1.0-lintx[k])*(1.0-linty[k]);
776: p2 = lintx[k] *(1.0-linty[k]);
777: p3 = lintx[k] * linty[k];
778: p4 = (1.0-lintx[k])* linty[k];
779: p = lintw[k]*PetscExpScalar( uLocal[0]*p1+ uLocal[1]*p2 + uLocal[2]*p3+ uLocal[3]*p4 );
780: int2 = int2 + p*p2;
781: }
782:
783: f[id] += rLocal[1] + int2*hxhy*(sc);
784:
785: /* element 3 */
786: uLocal[0] = x[j-1][i-1];
787: uLocal[1] = x[j-1][i];
788: uLocal[2] = x[j][i];
789: uLocal[3] = x[j][i-1];
790:
791: rLocal[2] = -1.0/3.0*uLocal[0]-1.0/6.0*uLocal[1]+2.0/3.0*uLocal[2]-1.0/6.0*uLocal[3];
792: rLocal[2] *= hxhy*alpha;
793: int3 = 0;
794: for(k=0; k<intnum; k++){
795: p1 = (1.0-lintx[k])*(1.0-linty[k]);
796: p2 = lintx[k] *(1.0-linty[k]);
797: p3 = lintx[k] * linty[k];
798: p4 = (1.0-lintx[k])* linty[k];
799: p = lintw[k]*PetscExpScalar( uLocal[0]*p1+ uLocal[1]*p2 + uLocal[2]*p3+ uLocal[3]*p4 );
800: int3 = int3 + p*p3;
801: }
802:
803: f[id] += rLocal[2] + int3*hxhy*(sc);
805: /* element 4 */
806: uLocal[0] = x[j-1][i];
807: uLocal[1] = x[j-1][i+1];
808: uLocal[2] = x[j][i+1];
809: uLocal[3] = x[j][i];
811: rLocal[3] = -1.0/6.0*uLocal[0]-1.0/3.0*uLocal[1]-1.0/6.0*uLocal[2]+2.0/3.0*uLocal[3];
812: rLocal[3] *= hxhy*alpha;
813: int4 = 0;
814: for(k=0; k<intnum; k++){
815: p1 = (1.0-lintx[k])*(1.0-linty[k]);
816: p2 = lintx[k] *(1.0-linty[k]);
817: p3 = lintx[k] * linty[k];
818: p4 = (1.0-lintx[k])* linty[k];
819: p = lintw[k]*PetscExpScalar( uLocal[0]*p1+ uLocal[1]*p2 + uLocal[2]*p3+ uLocal[3]*p4 );
820: int4 = int4 + p*p4;
821: }
822:
823: f[id] += rLocal[3] + int4*hxhy*(sc);
824: id++;
825: }
826: }
830: return(0);
831: }
833: //****************************************************************************