Actual source code: ex69.c
2: static char help[] = "Tests recovery from domain errors in MatMult() and PCApply()\n\n";
4: /*
5: See src/ksp/ksp/tutorials/ex19.c from which this was copied
6: */
8: #include <petscsnes.h>
9: #include <petscdm.h>
10: #include <petscdmda.h>
12: /*
13: User-defined routines and data structures
14: */
15: typedef struct {
16: PetscScalar u,v,omega,temp;
17: } Field;
19: PetscErrorCode FormFunctionLocal(DMDALocalInfo*,Field**,Field**,void*);
21: typedef struct {
22: PetscReal lidvelocity,prandtl,grashof; /* physical parameters */
23: PetscBool draw_contours; /* flag - 1 indicates drawing contours */
24: PetscBool errorindomain;
25: PetscBool errorindomainmf;
26: SNES snes;
27: } AppCtx;
29: typedef struct {
30: Mat Jmf;
31: } MatShellCtx;
33: extern PetscErrorCode FormInitialGuess(AppCtx*,DM,Vec);
34: extern PetscErrorCode MatMult_MyShell(Mat,Vec,Vec);
35: extern PetscErrorCode MatAssemblyEnd_MyShell(Mat,MatAssemblyType);
36: extern PetscErrorCode PCApply_MyShell(PC,Vec,Vec);
37: extern PetscErrorCode SNESComputeJacobian_MyShell(SNES,Vec,Mat,Mat,void*);
39: int main(int argc,char **argv)
40: {
41: AppCtx user; /* user-defined work context */
42: PetscInt mx,my;
44: MPI_Comm comm;
45: DM da;
46: Vec x;
47: Mat J = NULL,Jmf = NULL;
48: MatShellCtx matshellctx;
49: PetscInt mlocal,nlocal;
50: PC pc;
51: KSP ksp;
52: PetscBool errorinmatmult = PETSC_FALSE,errorinpcapply = PETSC_FALSE,errorinpcsetup = PETSC_FALSE;
54: PetscInitialize(&argc,&argv,(char*)0,help);
55: PetscOptionsGetBool(NULL,NULL,"-error_in_matmult",&errorinmatmult,NULL);
56: PetscOptionsGetBool(NULL,NULL,"-error_in_pcapply",&errorinpcapply,NULL);
57: PetscOptionsGetBool(NULL,NULL,"-error_in_pcsetup",&errorinpcsetup,NULL);
58: user.errorindomain = PETSC_FALSE;
59: PetscOptionsGetBool(NULL,NULL,"-error_in_domain",&user.errorindomain,NULL);
60: user.errorindomainmf = PETSC_FALSE;
61: PetscOptionsGetBool(NULL,NULL,"-error_in_domainmf",&user.errorindomainmf,NULL);
63: comm = PETSC_COMM_WORLD;
64: SNESCreate(comm,&user.snes);
66: /*
67: Create distributed array object to manage parallel grid and vectors
68: for principal unknowns (x) and governing residuals (f)
69: */
70: DMDACreate2d(PETSC_COMM_WORLD,DM_BOUNDARY_NONE,DM_BOUNDARY_NONE,DMDA_STENCIL_STAR,4,4,PETSC_DECIDE,PETSC_DECIDE,4,1,0,0,&da);
71: DMSetFromOptions(da);
72: DMSetUp(da);
73: SNESSetDM(user.snes,da);
75: DMDAGetInfo(da,0,&mx,&my,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,
76: PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE);
77: /*
78: Problem parameters (velocity of lid, prandtl, and grashof numbers)
79: */
80: user.lidvelocity = 1.0/(mx*my);
81: user.prandtl = 1.0;
82: user.grashof = 1.0;
84: PetscOptionsGetReal(NULL,NULL,"-lidvelocity",&user.lidvelocity,NULL);
85: PetscOptionsGetReal(NULL,NULL,"-prandtl",&user.prandtl,NULL);
86: PetscOptionsGetReal(NULL,NULL,"-grashof",&user.grashof,NULL);
87: PetscOptionsHasName(NULL,NULL,"-contours",&user.draw_contours);
89: DMDASetFieldName(da,0,"x_velocity");
90: DMDASetFieldName(da,1,"y_velocity");
91: DMDASetFieldName(da,2,"Omega");
92: DMDASetFieldName(da,3,"temperature");
94: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
95: Create user context, set problem data, create vector data structures.
96: Also, compute the initial guess.
97: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
99: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
100: Create nonlinear solver context
102: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
103: DMSetApplicationContext(da,&user);
104: DMDASNESSetFunctionLocal(da,INSERT_VALUES,(PetscErrorCode (*)(DMDALocalInfo*,void*,void*,void*))FormFunctionLocal,&user);
106: if (errorinmatmult) {
107: MatCreateSNESMF(user.snes,&Jmf);
108: MatSetFromOptions(Jmf);
109: MatGetLocalSize(Jmf,&mlocal,&nlocal);
110: matshellctx.Jmf = Jmf;
111: MatCreateShell(PetscObjectComm((PetscObject)Jmf),mlocal,nlocal,PETSC_DECIDE,PETSC_DECIDE,&matshellctx,&J);
112: MatShellSetOperation(J,MATOP_MULT,(void (*)(void))MatMult_MyShell);
113: MatShellSetOperation(J,MATOP_ASSEMBLY_END,(void (*)(void))MatAssemblyEnd_MyShell);
114: SNESSetJacobian(user.snes,J,J,MatMFFDComputeJacobian,NULL);
115: }
117: SNESSetFromOptions(user.snes);
118: PetscPrintf(comm,"lid velocity = %g, prandtl # = %g, grashof # = %g\n",(double)user.lidvelocity,(double)user.prandtl,(double)user.grashof);
120: if (errorinpcapply) {
121: SNESGetKSP(user.snes,&ksp);
122: KSPGetPC(ksp,&pc);
123: PCSetType(pc,PCSHELL);
124: PCShellSetApply(pc,PCApply_MyShell);
125: }
127: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
128: Solve the nonlinear system
129: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
130: DMCreateGlobalVector(da,&x);
131: FormInitialGuess(&user,da,x);
133: if (errorinpcsetup) {
134: SNESSetUp(user.snes);
135: SNESSetJacobian(user.snes,NULL,NULL,SNESComputeJacobian_MyShell,NULL);
136: }
137: SNESSolve(user.snes,NULL,x);
139: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
140: Free work space. All PETSc objects should be destroyed when they
141: are no longer needed.
142: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
143: MatDestroy(&J);
144: MatDestroy(&Jmf);
145: VecDestroy(&x);
146: DMDestroy(&da);
147: SNESDestroy(&user.snes);
148: PetscFinalize();
149: return 0;
150: }
152: /* ------------------------------------------------------------------- */
154: /*
155: FormInitialGuess - Forms initial approximation.
157: Input Parameters:
158: user - user-defined application context
159: X - vector
161: Output Parameter:
162: X - vector
163: */
164: PetscErrorCode FormInitialGuess(AppCtx *user,DM da,Vec X)
165: {
166: PetscInt i,j,mx,xs,ys,xm,ym;
167: PetscReal grashof,dx;
168: Field **x;
171: grashof = user->grashof;
173: DMDAGetInfo(da,0,&mx,0,0,0,0,0,0,0,0,0,0,0);
174: dx = 1.0/(mx-1);
176: /*
177: Get local grid boundaries (for 2-dimensional DMDA):
178: xs, ys - starting grid indices (no ghost points)
179: xm, ym - widths of local grid (no ghost points)
180: */
181: DMDAGetCorners(da,&xs,&ys,NULL,&xm,&ym,NULL);
183: /*
184: Get a pointer to vector data.
185: - For default PETSc vectors, VecGetArray() returns a pointer to
186: the data array. Otherwise, the routine is implementation dependent.
187: - You MUST call VecRestoreArray() when you no longer need access to
188: the array.
189: */
190: DMDAVecGetArray(da,X,&x);
192: /*
193: Compute initial guess over the locally owned part of the grid
194: Initial condition is motionless fluid and equilibrium temperature
195: */
196: for (j=ys; j<ys+ym; j++) {
197: for (i=xs; i<xs+xm; i++) {
198: x[j][i].u = 0.0;
199: x[j][i].v = 0.0;
200: x[j][i].omega = 0.0;
201: x[j][i].temp = (grashof>0)*i*dx;
202: }
203: }
205: /*
206: Restore vector
207: */
208: DMDAVecRestoreArray(da,X,&x);
209: return 0;
210: }
212: PetscErrorCode FormFunctionLocal(DMDALocalInfo *info,Field **x,Field **f,void *ptr)
213: {
214: AppCtx *user = (AppCtx*)ptr;
215: PetscInt xints,xinte,yints,yinte,i,j;
216: PetscReal hx,hy,dhx,dhy,hxdhy,hydhx;
217: PetscReal grashof,prandtl,lid;
218: PetscScalar u,uxx,uyy,vx,vy,avx,avy,vxp,vxm,vyp,vym;
219: static PetscInt fail = 0;
222: if ((fail++ > 7 && user->errorindomainmf) || (fail++ > 36 && user->errorindomain)) {
223: PetscMPIInt rank;
224: MPI_Comm_rank(PetscObjectComm((PetscObject)user->snes),&rank);
225: if (rank == 0) {
226: SNESSetFunctionDomainError(user->snes);
227: }
228: }
229: grashof = user->grashof;
230: prandtl = user->prandtl;
231: lid = user->lidvelocity;
233: /*
234: Define mesh intervals ratios for uniform grid.
236: Note: FD formulae below are normalized by multiplying through by
237: local volume element (i.e. hx*hy) to obtain coefficients O(1) in two dimensions.
239: */
240: dhx = (PetscReal)(info->mx-1); dhy = (PetscReal)(info->my-1);
241: hx = 1.0/dhx; hy = 1.0/dhy;
242: hxdhy = hx*dhy; hydhx = hy*dhx;
244: xints = info->xs; xinte = info->xs+info->xm; yints = info->ys; yinte = info->ys+info->ym;
246: /* Test whether we are on the bottom edge of the global array */
247: if (yints == 0) {
248: j = 0;
249: yints = yints + 1;
250: /* bottom edge */
251: for (i=info->xs; i<info->xs+info->xm; i++) {
252: f[j][i].u = x[j][i].u;
253: f[j][i].v = x[j][i].v;
254: f[j][i].omega = x[j][i].omega + (x[j+1][i].u - x[j][i].u)*dhy;
255: f[j][i].temp = x[j][i].temp-x[j+1][i].temp;
256: }
257: }
259: /* Test whether we are on the top edge of the global array */
260: if (yinte == info->my) {
261: j = info->my - 1;
262: yinte = yinte - 1;
263: /* top edge */
264: for (i=info->xs; i<info->xs+info->xm; i++) {
265: f[j][i].u = x[j][i].u - lid;
266: f[j][i].v = x[j][i].v;
267: f[j][i].omega = x[j][i].omega + (x[j][i].u - x[j-1][i].u)*dhy;
268: f[j][i].temp = x[j][i].temp-x[j-1][i].temp;
269: }
270: }
272: /* Test whether we are on the left edge of the global array */
273: if (xints == 0) {
274: i = 0;
275: xints = xints + 1;
276: /* left edge */
277: for (j=info->ys; j<info->ys+info->ym; j++) {
278: f[j][i].u = x[j][i].u;
279: f[j][i].v = x[j][i].v;
280: f[j][i].omega = x[j][i].omega - (x[j][i+1].v - x[j][i].v)*dhx;
281: f[j][i].temp = x[j][i].temp;
282: }
283: }
285: /* Test whether we are on the right edge of the global array */
286: if (xinte == info->mx) {
287: i = info->mx - 1;
288: xinte = xinte - 1;
289: /* right edge */
290: for (j=info->ys; j<info->ys+info->ym; j++) {
291: f[j][i].u = x[j][i].u;
292: f[j][i].v = x[j][i].v;
293: f[j][i].omega = x[j][i].omega - (x[j][i].v - x[j][i-1].v)*dhx;
294: f[j][i].temp = x[j][i].temp - (PetscReal)(grashof>0);
295: }
296: }
298: /* Compute over the interior points */
299: for (j=yints; j<yinte; j++) {
300: for (i=xints; i<xinte; i++) {
302: /*
303: convective coefficients for upwinding
304: */
305: vx = x[j][i].u; avx = PetscAbsScalar(vx);
306: vxp = .5*(vx+avx); vxm = .5*(vx-avx);
307: vy = x[j][i].v; avy = PetscAbsScalar(vy);
308: vyp = .5*(vy+avy); vym = .5*(vy-avy);
310: /* U velocity */
311: u = x[j][i].u;
312: uxx = (2.0*u - x[j][i-1].u - x[j][i+1].u)*hydhx;
313: uyy = (2.0*u - x[j-1][i].u - x[j+1][i].u)*hxdhy;
314: f[j][i].u = uxx + uyy - .5*(x[j+1][i].omega-x[j-1][i].omega)*hx;
316: /* V velocity */
317: u = x[j][i].v;
318: uxx = (2.0*u - x[j][i-1].v - x[j][i+1].v)*hydhx;
319: uyy = (2.0*u - x[j-1][i].v - x[j+1][i].v)*hxdhy;
320: f[j][i].v = uxx + uyy + .5*(x[j][i+1].omega-x[j][i-1].omega)*hy;
322: /* Omega */
323: u = x[j][i].omega;
324: uxx = (2.0*u - x[j][i-1].omega - x[j][i+1].omega)*hydhx;
325: uyy = (2.0*u - x[j-1][i].omega - x[j+1][i].omega)*hxdhy;
326: f[j][i].omega = uxx + uyy + (vxp*(u - x[j][i-1].omega) + vxm*(x[j][i+1].omega - u))*hy +
327: (vyp*(u - x[j-1][i].omega) + vym*(x[j+1][i].omega - u))*hx -
328: .5*grashof*(x[j][i+1].temp - x[j][i-1].temp)*hy;
330: /* Temperature */
331: u = x[j][i].temp;
332: uxx = (2.0*u - x[j][i-1].temp - x[j][i+1].temp)*hydhx;
333: uyy = (2.0*u - x[j-1][i].temp - x[j+1][i].temp)*hxdhy;
334: f[j][i].temp = uxx + uyy + prandtl*((vxp*(u - x[j][i-1].temp) + vxm*(x[j][i+1].temp - u))*hy +
335: (vyp*(u - x[j-1][i].temp) + vym*(x[j+1][i].temp - u))*hx);
336: }
337: }
339: /*
340: Flop count (multiply-adds are counted as 2 operations)
341: */
342: PetscLogFlops(84.0*info->ym*info->xm);
343: return 0;
344: }
346: PetscErrorCode MatMult_MyShell(Mat A,Vec x,Vec y)
347: {
348: MatShellCtx *matshellctx;
349: static PetscInt fail = 0;
351: MatShellGetContext(A,&matshellctx);
352: MatMult(matshellctx->Jmf,x,y);
353: if (fail++ > 5) {
354: PetscMPIInt rank;
355: MPI_Comm_rank(PetscObjectComm((PetscObject)A),&rank);
356: if (rank == 0) VecSetInf(y);
357: }
358: return 0;
359: }
361: PetscErrorCode MatAssemblyEnd_MyShell(Mat A,MatAssemblyType tp)
362: {
363: MatShellCtx *matshellctx;
365: MatShellGetContext(A,&matshellctx);
366: MatAssemblyEnd(matshellctx->Jmf,tp);
367: return 0;
368: }
370: PetscErrorCode PCApply_MyShell(PC pc,Vec x,Vec y)
371: {
372: static PetscInt fail = 0;
374: VecCopy(x,y);
375: if (fail++ > 3) {
376: PetscMPIInt rank;
377: MPI_Comm_rank(PetscObjectComm((PetscObject)pc),&rank);
378: if (rank == 0) VecSetInf(y);
379: }
380: return 0;
381: }
383: PETSC_EXTERN PetscErrorCode SNESComputeJacobian_DMDA(SNES,Vec,Mat,Mat,void*);
385: PetscErrorCode SNESComputeJacobian_MyShell(SNES snes,Vec X,Mat A,Mat B,void *ctx)
386: {
387: static PetscInt fail = 0;
389: SNESComputeJacobian_DMDA(snes,X,A,B,ctx);
390: if (fail++ > 0) {
391: MatZeroEntries(A);
392: }
393: return 0;
394: }
396: /*TEST
398: test:
399: args: -snes_converged_reason -ksp_converged_reason
401: test:
402: suffix: 2
403: args: -snes_converged_reason -ksp_converged_reason -error_in_matmult
405: test:
406: suffix: 3
407: args: -snes_converged_reason -ksp_converged_reason -error_in_pcapply
409: test:
410: suffix: 4
411: args: -snes_converged_reason -ksp_converged_reason -error_in_pcsetup
413: test:
414: suffix: 5
415: args: -snes_converged_reason -ksp_converged_reason -error_in_pcsetup -pc_type bjacobi
417: test:
418: suffix: 5_fieldsplit
419: args: -snes_converged_reason -ksp_converged_reason -error_in_pcsetup -pc_type fieldsplit
420: output_file: output/ex69_5.out
422: test:
423: suffix: 6
424: args: -snes_converged_reason -ksp_converged_reason -error_in_domainmf -snes_mf -pc_type none
426: test:
427: suffix: 7
428: args: -snes_converged_reason -ksp_converged_reason -error_in_domain
430: test:
431: suffix: 8
432: args: -snes_converged_reason -ksp_converged_reason -error_in_domain -snes_mf -pc_type none
434: TEST*/