Actual source code: ex5.c
petsc-3.7.7 2017-09-25
2: static char help[] = "Bratu nonlinear PDE in 2d.\n\
3: We solve the Bratu (SFI - solid fuel ignition) problem in a 2D rectangular\n\
4: domain, using distributed arrays (DMDAs) to partition the parallel grid.\n\
5: The command line options include:\n\
6: -par <parameter>, where <parameter> indicates the problem's nonlinearity\n\
7: problem SFI: <parameter> = Bratu parameter (0 <= par <= 6.81)\n\n\
8: -m_par/n_par <parameter>, where <parameter> indicates an integer\n \
9: that MMS3 will be evaluated with 2^m_par, 2^n_par";
11: /*T
12: Concepts: SNES^parallel Bratu example
13: Concepts: DMDA^using distributed arrays;
14: Concepts: IS coloirng types;
15: Processors: n
16: T*/
18: /* ------------------------------------------------------------------------
20: Solid Fuel Ignition (SFI) problem. This problem is modeled by
21: the partial differential equation
23: -Laplacian u - lambda*exp(u) = 0, 0 < x,y < 1,
25: with boundary conditions
27: u = 0 for x = 0, x = 1, y = 0, y = 1.
29: A finite difference approximation with the usual 5-point stencil
30: is used to discretize the boundary value problem to obtain a nonlinear
31: system of equations.
34: This example shows how geometric multigrid can be run transparently with a nonlinear solver so long
35: as SNESSetDM() is provided. Example usage
37: ./ex5 -pc_type mg -ksp_monitor -snes_view -pc_mg_levels 3 -pc_mg_galerkin -da_grid_x 17 -da_grid_y 17
38: -mg_levels_ksp_monitor -snes_monitor -mg_levels_pc_type sor -pc_mg_type full
40: or to run with grid sequencing on the nonlinear problem (note that you do not need to provide the number of
41: multigrid levels, it will be determined automatically based on the number of refinements done)
43: ./ex5 -pc_type mg -ksp_monitor -snes_view -pc_mg_galerkin -snes_grid_sequence 3
44: -mg_levels_ksp_monitor -snes_monitor -mg_levels_pc_type sor -pc_mg_type full
46: ------------------------------------------------------------------------- */
48: /*
49: Include "petscdmda.h" so that we can use distributed arrays (DMDAs).
50: Include "petscsnes.h" so that we can use SNES solvers. Note that this
51: */
52: #include <petscdm.h>
53: #include <petscdmda.h>
54: #include <petscsnes.h>
55: #include <petscmatlab.h>
57: /*
58: User-defined application context - contains data needed by the
59: application-provided call-back routines, FormJacobianLocal() and
60: FormFunctionLocal().
61: */
62: typedef struct {
63: PetscReal param; /* test problem parameter */
64: PetscInt mPar; /* MMS3 m parameter */
65: PetscInt nPar; /* MMS3 n parameter */
66: } AppCtx;
68: /*
69: User-defined routines
70: */
71: extern PetscErrorCode FormInitialGuess(DM,AppCtx*,Vec);
72: extern PetscErrorCode FormFunctionLocal(DMDALocalInfo*,PetscScalar**,PetscScalar**,AppCtx*);
73: extern PetscErrorCode FormExactSolution1(DM,AppCtx*,Vec);
74: extern PetscErrorCode FormFunctionLocalMMS1(DMDALocalInfo*,PetscScalar**,PetscScalar**,AppCtx*);
75: extern PetscErrorCode FormExactSolution2(DM,AppCtx*,Vec);
76: extern PetscErrorCode FormFunctionLocalMMS2(DMDALocalInfo*,PetscScalar**,PetscScalar**,AppCtx*);
77: extern PetscErrorCode FormExactSolution3(DM,AppCtx*,Vec);
78: extern PetscErrorCode FormFunctionLocalMMS3(DMDALocalInfo*,PetscScalar**,PetscScalar**,AppCtx*);
79: extern PetscErrorCode FormExactSolution4(DM,AppCtx*,Vec);
80: extern PetscErrorCode FormFunctionLocalMMS4(DMDALocalInfo*,PetscScalar**,PetscScalar**,AppCtx*);
81: extern PetscErrorCode FormJacobianLocal(DMDALocalInfo*,PetscScalar**,Mat,Mat,AppCtx*);
82: extern PetscErrorCode FormObjectiveLocal(DMDALocalInfo*,PetscScalar**,PetscReal*,AppCtx*);
83: #if defined(PETSC_HAVE_MATLAB_ENGINE)
84: extern PetscErrorCode FormFunctionMatlab(SNES,Vec,Vec,void*);
85: #endif
86: extern PetscErrorCode NonlinearGS(SNES,Vec,Vec,void*);
90: int main(int argc,char **argv)
91: {
92: SNES snes; /* nonlinear solver */
93: Vec x; /* solution vector */
94: AppCtx user; /* user-defined work context */
95: PetscInt its; /* iterations for convergence */
97: PetscReal bratu_lambda_max = 6.81;
98: PetscReal bratu_lambda_min = 0.;
99: PetscInt MMS = 0;
100: PetscBool flg = PETSC_FALSE;
101: DM da;
102: #if defined(PETSC_HAVE_MATLAB_ENGINE)
103: Vec r = NULL;
104: PetscBool matlab_function = PETSC_FALSE;
105: #endif
106: KSP ksp;
107: PetscInt lits,slits;
109: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
110: Initialize program
111: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
113: PetscInitialize(&argc,&argv,(char*)0,help);
115: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
116: Initialize problem parameters
117: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
118: user.param = 6.0;
119: PetscOptionsGetReal(NULL,NULL,"-par",&user.param,NULL);
120: if (user.param > bratu_lambda_max || user.param < bratu_lambda_min) SETERRQ3(PETSC_COMM_SELF,1,"Lambda, %g, is out of range, [%g, %g]", user.param, bratu_lambda_min, bratu_lambda_max);
121: PetscOptionsGetInt(NULL,NULL,"-mms",&MMS,NULL);
122: if (MMS == 3) {
123: user.mPar = 2;
124: user.nPar = 1;
125: PetscOptionsGetInt(NULL,NULL,"-m_par",&user.mPar,NULL);
126: PetscOptionsGetInt(NULL,NULL,"-n_par",&user.nPar,NULL);
127: }
129: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
130: Create nonlinear solver context
131: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
132: SNESCreate(PETSC_COMM_WORLD,&snes);
133: SNESSetCountersReset(snes,PETSC_FALSE);
134: SNESSetNGS(snes, NonlinearGS, NULL);
136: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
137: Create distributed array (DMDA) to manage parallel grid and vectors
138: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
139: DMDACreate2d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE,DMDA_STENCIL_STAR,-4,-4,PETSC_DECIDE,PETSC_DECIDE,1,1,NULL,NULL,&da);
140: DMDASetUniformCoordinates(da, 0.0, 1.0, 0.0, 1.0, 0.0, 1.0);
141: DMSetApplicationContext(da,&user);
142: SNESSetDM(snes,da);
143: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
144: Extract global vectors from DMDA; then duplicate for remaining
145: vectors that are the same types
146: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
147: DMCreateGlobalVector(da,&x);
149: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
150: Set local function evaluation routine
151: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
152: switch (MMS) {
153: case 1: DMDASNESSetFunctionLocal(da,INSERT_VALUES,(DMDASNESFunction)FormFunctionLocalMMS1,&user);break;
154: case 2: DMDASNESSetFunctionLocal(da,INSERT_VALUES,(DMDASNESFunction)FormFunctionLocalMMS2,&user);break;
155: case 3: DMDASNESSetFunctionLocal(da,INSERT_VALUES,(DMDASNESFunction)FormFunctionLocalMMS3,&user);break;
156: case 4: DMDASNESSetFunctionLocal(da,INSERT_VALUES,(DMDASNESFunction)FormFunctionLocalMMS4,&user);break;
157: default: DMDASNESSetFunctionLocal(da,INSERT_VALUES,(DMDASNESFunction)FormFunctionLocal,&user);
158: }
159: PetscOptionsGetBool(NULL,NULL,"-fd",&flg,NULL);
160: if (!flg) {
161: DMDASNESSetJacobianLocal(da,(DMDASNESJacobian)FormJacobianLocal,&user);
162: }
164: PetscOptionsGetBool(NULL,NULL,"-obj",&flg,NULL);
165: if (flg) {
166: DMDASNESSetObjectiveLocal(da,(DMDASNESObjective)FormObjectiveLocal,&user);
167: }
169: #if defined(PETSC_HAVE_MATLAB_ENGINE)
170: PetscOptionsGetBool(NULL,NULL,"-matlab_function",&matlab_function,0);
171: if (matlab_function) {
172: VecDuplicate(x,&r);
173: SNESSetFunction(snes,r,FormFunctionMatlab,&user);
174: }
175: #endif
177: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
178: Customize nonlinear solver; set runtime options
179: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
180: SNESSetFromOptions(snes);
182: FormInitialGuess(da,&user,x);
184: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
185: Solve nonlinear system
186: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
187: SNESSolve(snes,NULL,x);
188: SNESGetIterationNumber(snes,&its);
190: SNESGetLinearSolveIterations(snes,&slits);
191: SNESGetKSP(snes,&ksp);
192: KSPGetTotalIterations(ksp,&lits);
193: if (lits != slits) SETERRQ2(PETSC_COMM_WORLD,PETSC_ERR_PLIB,"Number of total linear iterations reported by SNES %D does not match reported by KSP %D",slits,lits);
194: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
195: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
197: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
198: If using MMS, check the l_2 error
199: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
200: if (MMS) {
201: Vec e;
202: PetscReal errorl2, errorinf;
203: PetscInt N;
205: VecDuplicate(x, &e);
206: PetscObjectViewFromOptions((PetscObject) x, NULL, "-sol_view");
207: if (MMS == 1) {FormExactSolution1(da, &user, e);}
208: else if (MMS == 2) {FormExactSolution2(da, &user, e);}
209: else if (MMS == 3) {FormExactSolution3(da, &user, e);}
210: else {FormExactSolution4(da, &user, e);}
211: PetscObjectViewFromOptions((PetscObject) e, NULL, "-exact_view");
212: VecAXPY(e, -1.0, x);
213: PetscObjectViewFromOptions((PetscObject) e, NULL, "-error_view");
214: VecNorm(e, NORM_2, &errorl2);
215: VecNorm(e, NORM_INFINITY, &errorinf);
216: VecGetSize(e, &N);
217: PetscPrintf(PETSC_COMM_WORLD, "N: %D error l2 %g inf %g\n", N, (double) errorl2/N, (double) errorinf);
218: VecDestroy(&e);
219: }
221: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
222: Free work space. All PETSc objects should be destroyed when they
223: are no longer needed.
224: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
225: #if defined(PETSC_HAVE_MATLAB_ENGINE)
226: VecDestroy(&r);
227: #endif
228: VecDestroy(&x);
229: SNESDestroy(&snes);
230: DMDestroy(&da);
231: PetscFinalize();
232: return 0;
233: }
234: /* ------------------------------------------------------------------- */
237: /*
238: FormInitialGuess - Forms initial approximation.
240: Input Parameters:
241: da - The DM
242: user - user-defined application context
244: Output Parameter:
245: X - vector
246: */
247: PetscErrorCode FormInitialGuess(DM da,AppCtx *user,Vec X)
248: {
249: PetscInt i,j,Mx,My,xs,ys,xm,ym;
251: PetscReal lambda,temp1,temp,hx,hy;
252: PetscScalar **x;
255: DMDAGetInfo(da,PETSC_IGNORE,&Mx,&My,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE);
257: lambda = user->param;
258: hx = 1.0/(PetscReal)(Mx-1);
259: hy = 1.0/(PetscReal)(My-1);
260: temp1 = lambda/(lambda + 1.0);
262: /*
263: Get a pointer to vector data.
264: - For default PETSc vectors, VecGetArray() returns a pointer to
265: the data array. Otherwise, the routine is implementation dependent.
266: - You MUST call VecRestoreArray() when you no longer need access to
267: the array.
268: */
269: DMDAVecGetArray(da,X,&x);
271: /*
272: Get local grid boundaries (for 2-dimensional DMDA):
273: xs, ys - starting grid indices (no ghost points)
274: xm, ym - widths of local grid (no ghost points)
276: */
277: DMDAGetCorners(da,&xs,&ys,NULL,&xm,&ym,NULL);
279: /*
280: Compute initial guess over the locally owned part of the grid
281: */
282: for (j=ys; j<ys+ym; j++) {
283: temp = (PetscReal)(PetscMin(j,My-j-1))*hy;
284: for (i=xs; i<xs+xm; i++) {
285: if (i == 0 || j == 0 || i == Mx-1 || j == My-1) {
286: /* boundary conditions are all zero Dirichlet */
287: x[j][i] = 0.0;
288: } else {
289: x[j][i] = temp1*PetscSqrtReal(PetscMin((PetscReal)(PetscMin(i,Mx-i-1))*hx,temp));
290: }
291: }
292: }
294: /*
295: Restore vector
296: */
297: DMDAVecRestoreArray(da,X,&x);
298: return(0);
299: }
303: /*
304: FormExactSolution1 - Forms initial approximation.
306: Input Parameters:
307: da - The DM
308: user - user-defined application context
310: Output Parameter:
311: X - vector
312: */
313: PetscErrorCode FormExactSolution1(DM da, AppCtx *user, Vec U)
314: {
315: DM coordDA;
316: Vec coordinates;
317: DMDACoor2d **coords;
318: PetscScalar **u;
319: PetscReal x, y;
320: PetscInt xs, ys, xm, ym, i, j;
324: DMDAGetCorners(da, &xs, &ys, NULL, &xm, &ym, NULL);
325: DMGetCoordinateDM(da, &coordDA);
326: DMGetCoordinates(da, &coordinates);
327: DMDAVecGetArray(coordDA, coordinates, &coords);
328: DMDAVecGetArray(da, U, &u);
329: for (j = ys; j < ys+ym; ++j) {
330: for (i = xs; i < xs+xm; ++i) {
331: x = PetscRealPart(coords[j][i].x);
332: y = PetscRealPart(coords[j][i].y);
333: u[j][i] = x*(1 - x)*y*(1 - y);
334: }
335: }
336: DMDAVecRestoreArray(da, U, &u);
337: DMDAVecRestoreArray(coordDA, coordinates, &coords);
338: return(0);
339: }
343: /*
344: FormExactSolution2 - Forms initial approximation.
346: Input Parameters:
347: da - The DM
348: user - user-defined application context
350: Output Parameter:
351: X - vector
352: */
353: PetscErrorCode FormExactSolution2(DM da, AppCtx *user, Vec U)
354: {
355: DM coordDA;
356: Vec coordinates;
357: DMDACoor2d **coords;
358: PetscScalar **u;
359: PetscReal x, y;
360: PetscInt xs, ys, xm, ym, i, j;
364: DMDAGetCorners(da, &xs, &ys, NULL, &xm, &ym, NULL);
365: DMGetCoordinateDM(da, &coordDA);
366: DMGetCoordinates(da, &coordinates);
367: DMDAVecGetArray(coordDA, coordinates, &coords);
368: DMDAVecGetArray(da, U, &u);
369: for (j = ys; j < ys+ym; ++j) {
370: for (i = xs; i < xs+xm; ++i) {
371: x = PetscRealPart(coords[j][i].x);
372: y = PetscRealPart(coords[j][i].y);
373: u[j][i] = PetscSinReal(PETSC_PI*x)*PetscSinReal(PETSC_PI*y);
374: }
375: }
376: DMDAVecRestoreArray(da, U, &u);
377: DMDAVecRestoreArray(coordDA, coordinates, &coords);
378: return(0);
379: }
383: /*
384: FormExactSolution3 - Forms initial approximation.
386: Input Parameters:
387: da - The DM
388: user - user-defined application context
390: Output Parameter:
391: X - vector
392: */
393: PetscErrorCode FormExactSolution3(DM da, AppCtx *user, Vec U)
394: {
395: DM coordDA;
396: Vec coordinates;
397: DMDACoor2d **coords;
398: PetscScalar **u;
399: PetscReal x, y;
400: PetscInt xs, ys, xm, ym, i, j, m, n;
403: m = PetscPowReal(2,user->mPar);
404: n = PetscPowReal(2,user->nPar);
407: DMDAGetCorners(da, &xs, &ys, NULL, &xm, &ym, NULL);
408: DMGetCoordinateDM(da, &coordDA);
409: DMGetCoordinates(da, &coordinates);
410: DMDAVecGetArray(coordDA, coordinates, &coords);
411: DMDAVecGetArray(da, U, &u);
413: for (j = ys; j < ys+ym; ++j) {
414: for (i = xs; i < xs+xm; ++i) {
415: x = PetscRealPart(coords[j][i].x);
416: y = PetscRealPart(coords[j][i].y);
418: u[j][i] = PetscSinReal(m*PETSC_PI*x*(1-y))*PetscSinReal(n*PETSC_PI*y*(1-x));
419: }
420: }
421: DMDAVecRestoreArray(da, U, &u);
422: DMDAVecRestoreArray(coordDA, coordinates, &coords);
423: return(0);
424: }
425: /* ------------------------------------------------------------------- */
429: /*
430: FormExactSolution4 - Forms initial approximation.
432: Input Parameters:
433: da - The DM
434: user - user-defined application context
436: Output Parameter:
437: X - vector
438: */
439: PetscErrorCode FormExactSolution4(DM da, AppCtx *user, Vec U)
440: {
441: DM coordDA;
442: Vec coordinates;
443: DMDACoor2d **coords;
444: PetscScalar **u;
445: PetscReal x, y, Lx, Ly;
446: PetscInt xs, ys, xm, ym, i, j;
450: DMDAGetCorners(da, &xs, &ys, NULL, &xm, &ym, NULL);
451: DMGetCoordinateDM(da, &coordDA);
452: DMGetCoordinates(da, &coordinates);
453: DMDAVecGetArray(coordDA, coordinates, &coords);
454: DMDAVecGetArray(da, U, &u);
456: Lx = PetscRealPart(coords[ys][xs+xm-1].x - coords[ys][xs].x);
457: Ly = PetscRealPart(coords[ys+ym-1][xs].y - coords[ys][xs].y);
459: for (j = ys; j < ys+ym; ++j) {
460: for (i = xs; i < xs+xm; ++i) {
461: x = PetscRealPart(coords[j][i].x);
462: y = PetscRealPart(coords[j][i].y);
463: u[j][i] = (PetscPowReal(x,4)-PetscSqr(Lx)*PetscSqr(x))*(PetscPowReal(y,4)-PetscSqr(Ly)*PetscSqr(y));
464: }
465: }
466: DMDAVecRestoreArray(da, U, &u);
467: DMDAVecRestoreArray(coordDA, coordinates, &coords);
468: return(0);
469: }
470: /* ------------------------------------------------------------------- */
473: /*
474: FormFunctionLocal - Evaluates nonlinear function, F(x) on local process patch
477: */
478: PetscErrorCode FormFunctionLocal(DMDALocalInfo *info,PetscScalar **x,PetscScalar **f,AppCtx *user)
479: {
481: PetscInt i,j;
482: PetscReal lambda,hx,hy,hxdhy,hydhx,sc;
483: PetscScalar u,ue,uw,un,us,uxx,uyy;
486: lambda = user->param;
487: hx = 1.0/(PetscReal)(info->mx-1);
488: hy = 1.0/(PetscReal)(info->my-1);
489: sc = hx*hy*lambda;
490: hxdhy = hx/hy;
491: hydhx = hy/hx;
492: /*
493: Compute function over the locally owned part of the grid
494: */
495: for (j=info->ys; j<info->ys+info->ym; j++) {
496: for (i=info->xs; i<info->xs+info->xm; i++) {
497: if (i == 0 || j == 0 || i == info->mx-1 || j == info->my-1) {
498: f[j][i] = 2.0*(hydhx+hxdhy)*x[j][i];
499: } else {
500: u = x[j][i];
501: uw = x[j][i-1];
502: ue = x[j][i+1];
503: un = x[j-1][i];
504: us = x[j+1][i];
506: if (i-1 == 0) uw = 0.;
507: if (i+1 == info->mx-1) ue = 0.;
508: if (j-1 == 0) un = 0.;
509: if (j+1 == info->my-1) us = 0.;
511: uxx = (2.0*u - uw - ue)*hydhx;
512: uyy = (2.0*u - un - us)*hxdhy;
513: f[j][i] = uxx + uyy - sc*PetscExpScalar(u);
514: }
515: }
516: }
517: PetscLogFlops(11.0*info->ym*info->xm);
518: return(0);
519: }
523: /* ---------------------------------------------------------------------------------
524: FormFunctionLocalMMS1 - Evaluates nonlinear function, F(x) on local process patch
526: u(x,y) = x(1-x)y(1-y)
528: -laplacian u* - lambda exp(u*) = 2x(1-x) + 2y(1-y) - lambda exp(x(1-x)y(1-y))
529:
530: Remark: the above is subtracted from the residual
531: -----------------------------------------------------------------------------------*/
532: PetscErrorCode FormFunctionLocalMMS1(DMDALocalInfo *info,PetscScalar **vx,PetscScalar **f,AppCtx *user)
533: {
535: PetscInt i,j;
536: PetscReal lambda,hx,hy,hxdhy,hydhx;
537: PetscScalar u,ue,uw,un,us,uxx,uyy;
538: PetscReal x,y;
539: DM coordDA;
540: Vec coordinates;
541: DMDACoor2d **coords;
542: Vec bcv = NULL;
543: PetscScalar **bcx = NULL;
546: lambda = user->param;
547: /* Extract coordinates */
548: DMGetCoordinateDM(info->da, &coordDA);
549: DMGetCoordinates(info->da, &coordinates);
550: DMDAVecGetArray(coordDA, coordinates, &coords);
551: hx = info->xm > 1 ? PetscRealPart(coords[info->ys][info->xs+1].x) - PetscRealPart(coords[info->ys][info->xs].x) : 1.0;
552: hy = info->ym > 1 ? PetscRealPart(coords[info->ys+1][info->xs].y) - PetscRealPart(coords[info->ys][info->xs].y) : 1.0;
553: hxdhy = hx/hy;
554: hydhx = hy/hx;
555: DMGetNamedLocalVector(info->da, "_petsc_boundary_conditions_", &bcv);
556: DMDAVecGetArray(info->da, bcv, &bcx);
557: /* Compute function over the locally owned part of the grid */
558: for (j=info->ys; j<info->ys+info->ym; j++) {
559: for (i=info->xs; i<info->xs+info->xm; i++) {
560: if (i == 0 || j == 0 || i == info->mx-1 || j == info->my-1) {
561: f[j][i] = 2.0*(hydhx+hxdhy)*(vx[j][i] - bcx[j][i]);
562: } else {
563: x = PetscRealPart(coords[j][i].x);
564: y = PetscRealPart(coords[j][i].y);
565: u = vx[j][i];
566: uw = vx[j][i-1];
567: ue = vx[j][i+1];
568: un = vx[j-1][i];
569: us = vx[j+1][i];
571: if (i-1 == 0) uw = bcx[j][i-1];
572: if (i+1 == info->mx-1) ue = bcx[j][i+1];
573: if (j-1 == 0) un = bcx[j-1][i];
574: if (j+1 == info->my-1) us = bcx[j+1][i];
576: uxx = (2.0*u - uw - ue)*hydhx;
577: uyy = (2.0*u - un - us)*hxdhy;
578: f[j][i] = uxx + uyy - hx*hy*(lambda*PetscExpScalar(u) + 2*x*(1 - x) + 2*y*(1 - y) - lambda*exp(x*(1 - x)*y*(1 - y)));
579: }
580: }
581: }
582: DMDAVecRestoreArray(info->da, bcv, &bcx);
583: DMRestoreNamedLocalVector(info->da, "_petsc_boundary_conditions_", &bcv);
584: DMDAVecRestoreArray(coordDA, coordinates, &coords);
585: PetscLogFlops(11.0*info->ym*info->xm);
586: return(0);
587: }
591: /* ---------------------------------------------------------------------------------
592: FormFunctionLocalMMS2 - Evaluates nonlinear function, F(x) on local process patch
594: u(x,y) = sin(pi x)sin(pi y)
596: -laplacian u* - lambda exp(u*) = 2 pi^2 sin(pi x) sin(pi y) - lambda exp(sin(pi x)sin(pi y))
598: Remark: the above is subtracted from the residual
599: -----------------------------------------------------------------------------------*/
600: PetscErrorCode FormFunctionLocalMMS2(DMDALocalInfo *info,PetscScalar **vx,PetscScalar **f,AppCtx *user)
601: {
603: PetscInt i,j;
604: PetscReal lambda,hx,hy,hxdhy,hydhx;
605: PetscScalar u,ue,uw,un,us,uxx,uyy;
606: PetscReal x,y;
607: DM coordDA;
608: Vec coordinates;
609: DMDACoor2d **coords;
610: Vec bcv = NULL;
611: PetscScalar **bcx = NULL;
614: lambda = user->param;
615: /* Extract coordinates */
616: DMGetCoordinateDM(info->da, &coordDA);
617: DMGetCoordinates(info->da, &coordinates);
618: DMDAVecGetArray(coordDA, coordinates, &coords);
619: hx = info->xm > 1 ? PetscRealPart(coords[info->ys][info->xs+1].x) - PetscRealPart(coords[info->ys][info->xs].x) : 1.0;
620: hy = info->ym > 1 ? PetscRealPart(coords[info->ys+1][info->xs].y) - PetscRealPart(coords[info->ys][info->xs].y) : 1.0;
621: hxdhy = hx/hy;
622: hydhx = hy/hx;
623: DMGetNamedLocalVector(info->da, "_petsc_boundary_conditions_", &bcv);
624: DMDAVecGetArray(info->da, bcv, &bcx);
625: /* Compute function over the locally owned part of the grid */
626: for (j=info->ys; j<info->ys+info->ym; j++) {
627: for (i=info->xs; i<info->xs+info->xm; i++) {
628: if (i == 0 || j == 0 || i == info->mx-1 || j == info->my-1) {
629: f[j][i] = 2.0*(hydhx+hxdhy)*(vx[j][i] - bcx[j][i]);
630: } else {
631: x = PetscRealPart(coords[j][i].x);
632: y = PetscRealPart(coords[j][i].y);
633: u = vx[j][i];
634: uw = vx[j][i-1];
635: ue = vx[j][i+1];
636: un = vx[j-1][i];
637: us = vx[j+1][i];
639: if (i-1 == 0) uw = bcx[j][i-1];
640: if (i+1 == info->mx-1) ue = bcx[j][i+1];
641: if (j-1 == 0) un = bcx[j-1][i];
642: if (j+1 == info->my-1) us = bcx[j+1][i];
644: uxx = (2.0*u - uw - ue)*hydhx;
645: uyy = (2.0*u - un - us)*hxdhy;
646: f[j][i] = uxx + uyy - hx*hy*(lambda*PetscExpScalar(u) + 2*PetscSqr(PETSC_PI)*PetscSinReal(PETSC_PI*x)*PetscSinReal(PETSC_PI*y) - lambda*exp(PetscSinReal(PETSC_PI*x)*PetscSinReal(PETSC_PI*y)));
647: }
648: }
649: }
650: DMDAVecRestoreArray(info->da, bcv, &bcx);
651: DMRestoreNamedLocalVector(info->da, "_petsc_boundary_conditions_", &bcv);
652: DMDAVecRestoreArray(coordDA, coordinates, &coords);
653: PetscLogFlops(11.0*info->ym*info->xm);
654: return(0);
655: }
659: /* ---------------------------------------------------------------------------------
660: FormFunctionLocalMMS3 - Evaluates nonlinear function, F(x) on local process patch
662: u(x,y) = sin(m pi x(1-y))sin(n pi y(1-x))
664: -laplacian u* - lambda exp(u*) = -(exp(Sin[m*Pi*x*(1 - y)]*Sin[n*Pi*(1 - x)*y])*lambda) +
665: Pi^2*(-2*m*n*((-1 + x)*x + (-1 + y)*y)*Cos[m*Pi*x*(-1 + y)]*
666: Cos[n*Pi*(-1 + x)*y] + (m^2*(x^2 + (-1 + y)^2) + n^2*((-1 + x)^2 + y^2))*
667: Sin[m*Pi*x*(-1 + y)]*Sin[n*Pi*(-1 + x)*y])
669: Remark: the above is subtracted from the residual
670: -----------------------------------------------------------------------------------*/
671: PetscErrorCode FormFunctionLocalMMS3(DMDALocalInfo *info,PetscScalar **vx,PetscScalar **f,AppCtx *user)
672: {
674: PetscInt i,j;
675: PetscReal lambda,hx,hy,hxdhy,hydhx,m,n;
676: PetscScalar u,ue,uw,un,us,uxx,uyy;
677: PetscReal x,y;
678: DM coordDA;
679: Vec coordinates;
680: DMDACoor2d **coords;
682: m = PetscPowReal(2,user->mPar);
683: n = PetscPowReal(2,user->nPar);
686: lambda = user->param;
687: hx = 1.0/(PetscReal)(info->mx-1);
688: hy = 1.0/(PetscReal)(info->my-1);
689: hxdhy = hx/hy;
690: hydhx = hy/hx;
691: /* Extract coordinates */
692: DMGetCoordinateDM(info->da, &coordDA);
693: DMGetCoordinates(info->da, &coordinates);
694: DMDAVecGetArray(coordDA, coordinates, &coords);
695: /* Compute function over the locally owned part of the grid */
696: for (j=info->ys; j<info->ys+info->ym; j++) {
697: for (i=info->xs; i<info->xs+info->xm; i++) {
698: if (i == 0 || j == 0 || i == info->mx-1 || j == info->my-1) {
699: f[j][i] = 2.0*(hydhx+hxdhy)*vx[j][i];
700: } else {
701: x = PetscRealPart(coords[j][i].x);
702: y = PetscRealPart(coords[j][i].y);
703: u = vx[j][i];
704: uw = vx[j][i-1];
705: ue = vx[j][i+1];
706: un = vx[j-1][i];
707: us = vx[j+1][i];
709: if (i-1 == 0) uw = 0.;
710: if (i+1 == info->mx-1) ue = 0.;
711: if (j-1 == 0) un = 0.;
712: if (j+1 == info->my-1) us = 0.;
714: uxx = (2.0*u - uw - ue)*hydhx;
715: uyy = (2.0*u - un - us)*hxdhy;
717: f[j][i] = uxx + uyy - hx*hy*(lambda*PetscExpScalar(u)-(PetscExpReal(PetscSinReal(m*PETSC_PI*x*(1 - y))*PetscSinReal(n*PETSC_PI*(1 - x)*y))*lambda) + PetscSqr(PETSC_PI)*(-2*m*n*((-1 + x)*x + (-1 + y)*y)*PetscCosReal(m*PETSC_PI*x*(-1 + y))*PetscCosReal(n*PETSC_PI*(-1 + x)*y) + (PetscSqr(m)*(PetscSqr(x) + PetscSqr(-1 + y)) + PetscSqr(n)*(PetscSqr(-1 + x) + PetscSqr(y)))*PetscSinReal(m*PETSC_PI*x*(-1 + y))*PetscSinReal(n*PETSC_PI*(-1 + x)*y)));
718: }
719: }
720: }
721: DMDAVecRestoreArray(coordDA, coordinates, &coords);
722: PetscLogFlops(11.0*info->ym*info->xm);
723: return(0);
725: }
729: /* ---------------------------------------------------------------------------------
730: FormFunctionLocalMMS4 - Evaluates nonlinear function, F(x) on local process patch
732: u(x,y) = (x^4 - Lx^2 x^2)(y^2-Ly^2 y^2)
734: -laplacian u* - lambda exp(u*) = 2x^2(x^2-Lx^2)(Ly^2-6y^2)+2y^2(Lx^2-6x^2)(y^2-Ly^2)
735: -lambda exp((x^4 - Lx^2 x^2)(y^2-Ly^2 y^2))
737: Remark: the above is subtracted from the residual
738: -----------------------------------------------------------------------------------*/
739: PetscErrorCode FormFunctionLocalMMS4(DMDALocalInfo *info,PetscScalar **vx,PetscScalar **f,AppCtx *user)
740: {
742: PetscInt i,j;
743: PetscReal lambda,hx,hy,hxdhy,hydhx,Lx,Ly;
744: PetscScalar u,ue,uw,un,us,uxx,uyy;
745: PetscReal x,y;
746: DM coordDA;
747: Vec coordinates;
748: DMDACoor2d **coords;
751: lambda = user->param;
752: hx = 1.0/(PetscReal)(info->mx-1);
753: hy = 1.0/(PetscReal)(info->my-1);
754: hxdhy = hx/hy;
755: hydhx = hy/hx;
756: /* Extract coordinates */
757: DMGetCoordinateDM(info->da, &coordDA);
758: DMGetCoordinates(info->da, &coordinates);
759: DMDAVecGetArray(coordDA, coordinates, &coords);
761: Lx = PetscRealPart(coords[info->ys][info->xs+info->xm-1].x - coords[info->ys][info->xs].x);
762: Ly = PetscRealPart(coords[info->ys+info->ym-1][info->xs].y - coords[info->ys][info->xs].y);
764: /* Compute function over the locally owned part of the grid */
765: for (j=info->ys; j<info->ys+info->ym; j++) {
766: for (i=info->xs; i<info->xs+info->xm; i++) {
767: if (i == 0 || j == 0 || i == info->mx-1 || j == info->my-1) {
768: f[j][i] = 2.0*(hydhx+hxdhy)*vx[j][i];
769: } else {
770: x = PetscRealPart(coords[j][i].x);
771: y = PetscRealPart(coords[j][i].y);
772: u = vx[j][i];
773: uw = vx[j][i-1];
774: ue = vx[j][i+1];
775: un = vx[j-1][i];
776: us = vx[j+1][i];
778: if (i-1 == 0) uw = 0.;
779: if (i+1 == info->mx-1) ue = 0.;
780: if (j-1 == 0) un = 0.;
781: if (j+1 == info->my-1) us = 0.;
783: uxx = (2.0*u - uw - ue)*hydhx;
784: uyy = (2.0*u - un - us)*hxdhy;
785: f[j][i] = uxx + uyy - hx*hy*(lambda*PetscExpScalar(u)+ 2*PetscSqr(x)*(PetscSqr(x)-PetscSqr(Lx))*(PetscSqr(Ly)-6*PetscSqr(y))+2*PetscSqr(y)*(PetscSqr(Lx)-6*PetscSqr(x))*(PetscSqr(y)-PetscSqr(Ly))-lambda*PetscExpReal((PetscPowReal(x,4)-PetscSqr(Lx)*PetscSqr(x))*(PetscPowReal(y,4)-PetscSqr(Ly)*PetscSqr(y))));
786: }
787: }
788: }
789: DMDAVecRestoreArray(coordDA, coordinates, &coords);
790: PetscLogFlops(11.0*info->ym*info->xm);
791: return(0);
793: }
797: /* FormObjectiveLocal - Evaluates nonlinear function, F(x) on local process patch */
798: PetscErrorCode FormObjectiveLocal(DMDALocalInfo *info,PetscScalar **x,PetscReal *obj,AppCtx *user)
799: {
801: PetscInt i,j;
802: PetscReal lambda,hx,hy,hxdhy,hydhx,sc,lobj=0;
803: PetscScalar u,ue,uw,un,us,uxux,uyuy;
804: MPI_Comm comm;
807: *obj = 0;
808: PetscObjectGetComm((PetscObject)info->da,&comm);
809: lambda = user->param;
810: hx = 1.0/(PetscReal)(info->mx-1);
811: hy = 1.0/(PetscReal)(info->my-1);
812: sc = hx*hy*lambda;
813: hxdhy = hx/hy;
814: hydhx = hy/hx;
815: /*
816: Compute function over the locally owned part of the grid
817: */
818: for (j=info->ys; j<info->ys+info->ym; j++) {
819: for (i=info->xs; i<info->xs+info->xm; i++) {
820: if (i == 0 || j == 0 || i == info->mx-1 || j == info->my-1) {
821: lobj += PetscRealPart((hydhx + hxdhy)*x[j][i]*x[j][i]);
822: } else {
823: u = x[j][i];
824: uw = x[j][i-1];
825: ue = x[j][i+1];
826: un = x[j-1][i];
827: us = x[j+1][i];
829: if (i-1 == 0) uw = 0.;
830: if (i+1 == info->mx-1) ue = 0.;
831: if (j-1 == 0) un = 0.;
832: if (j+1 == info->my-1) us = 0.;
834: /* F[u] = 1/2\int_{\omega}\nabla^2u(x)*u(x)*dx */
836: uxux = u*(2.*u - ue - uw)*hydhx;
837: uyuy = u*(2.*u - un - us)*hxdhy;
839: lobj += PetscRealPart(0.5*(uxux + uyuy) - sc*PetscExpScalar(u));
840: }
841: }
842: }
843: PetscLogFlops(12.0*info->ym*info->xm);
844: MPI_Allreduce(&lobj,obj,1,MPIU_REAL,MPIU_SUM,comm);
845: return(0);
846: }
850: /*
851: FormJacobianLocal - Evaluates Jacobian matrix on local process patch
852: */
853: PetscErrorCode FormJacobianLocal(DMDALocalInfo *info,PetscScalar **x,Mat jac,Mat jacpre,AppCtx *user)
854: {
856: PetscInt i,j,k;
857: MatStencil col[5],row;
858: PetscScalar lambda,v[5],hx,hy,hxdhy,hydhx,sc;
859: DM coordDA;
860: Vec coordinates;
861: DMDACoor2d **coords;
864: lambda = user->param;
865: /* Extract coordinates */
866: DMGetCoordinateDM(info->da, &coordDA);
867: DMGetCoordinates(info->da, &coordinates);
868: DMDAVecGetArray(coordDA, coordinates, &coords);
869: hx = info->xm > 1 ? PetscRealPart(coords[info->ys][info->xs+1].x) - PetscRealPart(coords[info->ys][info->xs].x) : 1.0;
870: hy = info->ym > 1 ? PetscRealPart(coords[info->ys+1][info->xs].y) - PetscRealPart(coords[info->ys][info->xs].y) : 1.0;
871: DMDAVecRestoreArray(coordDA, coordinates, &coords);
872: hxdhy = hx/hy;
873: hydhx = hy/hx;
874: sc = hx*hy*lambda;
877: /*
878: Compute entries for the locally owned part of the Jacobian.
879: - Currently, all PETSc parallel matrix formats are partitioned by
880: contiguous chunks of rows across the processors.
881: - Each processor needs to insert only elements that it owns
882: locally (but any non-local elements will be sent to the
883: appropriate processor during matrix assembly).
884: - Here, we set all entries for a particular row at once.
885: - We can set matrix entries either using either
886: MatSetValuesLocal() or MatSetValues(), as discussed above.
887: */
888: for (j=info->ys; j<info->ys+info->ym; j++) {
889: for (i=info->xs; i<info->xs+info->xm; i++) {
890: row.j = j; row.i = i;
891: /* boundary points */
892: if (i == 0 || j == 0 || i == info->mx-1 || j == info->my-1) {
893: v[0] = 2.0*(hydhx + hxdhy);
894: MatSetValuesStencil(jacpre,1,&row,1,&row,v,INSERT_VALUES);
895: } else {
896: k = 0;
897: /* interior grid points */
898: if (j-1 != 0) {
899: v[k] = -hxdhy;
900: col[k].j = j - 1; col[k].i = i;
901: k++;
902: }
903: if (i-1 != 0) {
904: v[k] = -hydhx;
905: col[k].j = j; col[k].i = i-1;
906: k++;
907: }
909: v[k] = 2.0*(hydhx + hxdhy) - sc*PetscExpScalar(x[j][i]); col[k].j = row.j; col[k].i = row.i; k++;
911: if (i+1 != info->mx-1) {
912: v[k] = -hydhx;
913: col[k].j = j; col[k].i = i+1;
914: k++;
915: }
916: if (j+1 != info->mx-1) {
917: v[k] = -hxdhy;
918: col[k].j = j + 1; col[k].i = i;
919: k++;
920: }
921: MatSetValuesStencil(jacpre,1,&row,k,col,v,INSERT_VALUES);
922: }
923: }
924: }
926: /*
927: Assemble matrix, using the 2-step process:
928: MatAssemblyBegin(), MatAssemblyEnd().
929: */
930: MatAssemblyBegin(jacpre,MAT_FINAL_ASSEMBLY);
931: MatAssemblyEnd(jacpre,MAT_FINAL_ASSEMBLY);
932: /*
933: Tell the matrix we will never add a new nonzero location to the
934: matrix. If we do, it will generate an error.
935: */
936: MatSetOption(jac,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);
937: return(0);
938: }
940: #if defined(PETSC_HAVE_MATLAB_ENGINE)
943: PetscErrorCode FormFunctionMatlab(SNES snes,Vec X,Vec F,void *ptr)
944: {
945: AppCtx *user = (AppCtx*)ptr;
947: PetscInt Mx,My;
948: PetscReal lambda,hx,hy;
949: Vec localX,localF;
950: MPI_Comm comm;
951: DM da;
954: SNESGetDM(snes,&da);
955: DMGetLocalVector(da,&localX);
956: DMGetLocalVector(da,&localF);
957: PetscObjectSetName((PetscObject)localX,"localX");
958: PetscObjectSetName((PetscObject)localF,"localF");
959: DMDAGetInfo(da,PETSC_IGNORE,&Mx,&My,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE);
961: lambda = user->param;
962: hx = 1.0/(PetscReal)(Mx-1);
963: hy = 1.0/(PetscReal)(My-1);
965: PetscObjectGetComm((PetscObject)snes,&comm);
966: /*
967: Scatter ghost points to local vector,using the 2-step process
968: DMGlobalToLocalBegin(),DMGlobalToLocalEnd().
969: By placing code between these two statements, computations can be
970: done while messages are in transition.
971: */
972: DMGlobalToLocalBegin(da,X,INSERT_VALUES,localX);
973: DMGlobalToLocalEnd(da,X,INSERT_VALUES,localX);
974: PetscMatlabEnginePut(PETSC_MATLAB_ENGINE_(comm),(PetscObject)localX);
975: PetscMatlabEngineEvaluate(PETSC_MATLAB_ENGINE_(comm),"localF=ex5m(localX,%18.16e,%18.16e,%18.16e)",hx,hy,lambda);
976: PetscMatlabEngineGet(PETSC_MATLAB_ENGINE_(comm),(PetscObject)localF);
978: /*
979: Insert values into global vector
980: */
981: DMLocalToGlobalBegin(da,localF,INSERT_VALUES,F);
982: DMLocalToGlobalEnd(da,localF,INSERT_VALUES,F);
983: DMRestoreLocalVector(da,&localX);
984: DMRestoreLocalVector(da,&localF);
985: return(0);
986: }
987: #endif
989: /* ------------------------------------------------------------------- */
992: /*
993: Applies some sweeps on nonlinear Gauss-Seidel on each process
995: */
996: PetscErrorCode NonlinearGS(SNES snes,Vec X, Vec B, void *ctx)
997: {
998: PetscInt i,j,k,Mx,My,xs,ys,xm,ym,its,tot_its,sweeps,l;
1000: PetscReal lambda,hx,hy,hxdhy,hydhx,sc;
1001: PetscScalar **x,**b,bij,F,F0=0,J,u,un,us,ue,eu,uw,uxx,uyy,y;
1002: PetscReal atol,rtol,stol;
1003: DM da;
1004: AppCtx *user;
1005: Vec localX,localB;
1008: tot_its = 0;
1009: SNESNGSGetSweeps(snes,&sweeps);
1010: SNESNGSGetTolerances(snes,&atol,&rtol,&stol,&its);
1011: SNESGetDM(snes,&da);
1012: DMGetApplicationContext(da,(void**)&user);
1014: DMDAGetInfo(da,PETSC_IGNORE,&Mx,&My,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE);
1016: lambda = user->param;
1017: hx = 1.0/(PetscReal)(Mx-1);
1018: hy = 1.0/(PetscReal)(My-1);
1019: sc = hx*hy*lambda;
1020: hxdhy = hx/hy;
1021: hydhx = hy/hx;
1024: DMGetLocalVector(da,&localX);
1025: if (B) {
1026: DMGetLocalVector(da,&localB);
1027: }
1028: for (l=0; l<sweeps; l++) {
1030: DMGlobalToLocalBegin(da,X,INSERT_VALUES,localX);
1031: DMGlobalToLocalEnd(da,X,INSERT_VALUES,localX);
1032: if (B) {
1033: DMGlobalToLocalBegin(da,B,INSERT_VALUES,localB);
1034: DMGlobalToLocalEnd(da,B,INSERT_VALUES,localB);
1035: }
1036: /*
1037: Get a pointer to vector data.
1038: - For default PETSc vectors, VecGetArray() returns a pointer to
1039: the data array. Otherwise, the routine is implementation dependent.
1040: - You MUST call VecRestoreArray() when you no longer need access to
1041: the array.
1042: */
1043: DMDAVecGetArray(da,localX,&x);
1044: if (B) DMDAVecGetArray(da,localB,&b);
1045: /*
1046: Get local grid boundaries (for 2-dimensional DMDA):
1047: xs, ys - starting grid indices (no ghost points)
1048: xm, ym - widths of local grid (no ghost points)
1049: */
1050: DMDAGetCorners(da,&xs,&ys,NULL,&xm,&ym,NULL);
1052: for (j=ys; j<ys+ym; j++) {
1053: for (i=xs; i<xs+xm; i++) {
1054: if (i == 0 || j == 0 || i == Mx-1 || j == My-1) {
1055: /* boundary conditions are all zero Dirichlet */
1056: x[j][i] = 0.0;
1057: } else {
1058: if (B) bij = b[j][i];
1059: else bij = 0.;
1061: u = x[j][i];
1062: un = x[j-1][i];
1063: us = x[j+1][i];
1064: ue = x[j][i-1];
1065: uw = x[j][i+1];
1067: for (k=0; k<its; k++) {
1068: eu = PetscExpScalar(u);
1069: uxx = (2.0*u - ue - uw)*hydhx;
1070: uyy = (2.0*u - un - us)*hxdhy;
1071: F = uxx + uyy - sc*eu - bij;
1072: if (k == 0) F0 = F;
1073: J = 2.0*(hydhx + hxdhy) - sc*eu;
1074: y = F/J;
1075: u -= y;
1076: tot_its++;
1078: if (atol > PetscAbsReal(PetscRealPart(F)) ||
1079: rtol*PetscAbsReal(PetscRealPart(F0)) > PetscAbsReal(PetscRealPart(F)) ||
1080: stol*PetscAbsReal(PetscRealPart(u)) > PetscAbsReal(PetscRealPart(y))) {
1081: break;
1082: }
1083: }
1084: x[j][i] = u;
1085: }
1086: }
1087: }
1088: /*
1089: Restore vector
1090: */
1091: DMDAVecRestoreArray(da,localX,&x);
1092: DMLocalToGlobalBegin(da,localX,INSERT_VALUES,X);
1093: DMLocalToGlobalEnd(da,localX,INSERT_VALUES,X);
1094: }
1095: PetscLogFlops(tot_its*(21.0));
1096: DMRestoreLocalVector(da,&localX);
1097: if (B) {
1098: DMDAVecRestoreArray(da,localB,&b);
1099: DMRestoreLocalVector(da,&localB);
1100: }
1101: return(0);
1102: }