Actual source code: ex69.c
petsc-3.7.3 2016-08-01
2: static char help[] = "Tests recovery from domain errors in MatMult() and PCApply()\n\n";
4: /*
5: See src/ksp/ksp/examples/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: PETSC_EXTERN PetscErrorCode SNESComputeJacobian_MyShell(SNES,Vec,Mat,Mat,void*);
41: int main(int argc,char **argv)
42: {
43: AppCtx user; /* user-defined work context */
44: PetscInt mx,my;
46: MPI_Comm comm;
47: DM da;
48: Vec x;
49: Mat J = NULL,Jmf = NULL;
50: MatShellCtx matshellctx;
51: PetscInt mlocal,nlocal;
52: PC pc;
53: KSP ksp;
54: PetscBool errorinmatmult = PETSC_FALSE,errorinpcapply = PETSC_FALSE,errorinpcsetup = PETSC_FALSE;
56: PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return(1);
59: PetscOptionsGetBool(NULL,NULL,"-error_in_matmult",&errorinmatmult,NULL);
60: PetscOptionsGetBool(NULL,NULL,"-error_in_pcapply",&errorinpcapply,NULL);
61: PetscOptionsGetBool(NULL,NULL,"-error_in_pcsetup",&errorinpcsetup,NULL);
62: user.errorindomain = PETSC_FALSE;
63: PetscOptionsGetBool(NULL,NULL,"-error_in_domain",&user.errorindomain,NULL);
64: user.errorindomainmf = PETSC_FALSE;
65: PetscOptionsGetBool(NULL,NULL,"-error_in_domainmf",&user.errorindomainmf,NULL);
67: comm = PETSC_COMM_WORLD;
68: SNESCreate(comm,&user.snes);
70: /*
71: Create distributed array object to manage parallel grid and vectors
72: for principal unknowns (x) and governing residuals (f)
73: */
74: DMDACreate2d(PETSC_COMM_WORLD,DM_BOUNDARY_NONE,DM_BOUNDARY_NONE,DMDA_STENCIL_STAR,-4,-4,PETSC_DECIDE,PETSC_DECIDE,4,1,0,0,&da);
75: SNESSetDM(user.snes,da);
77: DMDAGetInfo(da,0,&mx,&my,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,
78: PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE);
79: /*
80: Problem parameters (velocity of lid, prandtl, and grashof numbers)
81: */
82: user.lidvelocity = 1.0/(mx*my);
83: user.prandtl = 1.0;
84: user.grashof = 1.0;
86: PetscOptionsGetReal(NULL,NULL,"-lidvelocity",&user.lidvelocity,NULL);
87: PetscOptionsGetReal(NULL,NULL,"-prandtl",&user.prandtl,NULL);
88: PetscOptionsGetReal(NULL,NULL,"-grashof",&user.grashof,NULL);
89: PetscOptionsHasName(NULL,NULL,"-contours",&user.draw_contours);
91: DMDASetFieldName(da,0,"x_velocity");
92: DMDASetFieldName(da,1,"y_velocity");
93: DMDASetFieldName(da,2,"Omega");
94: DMDASetFieldName(da,3,"temperature");
96: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
97: Create user context, set problem data, create vector data structures.
98: Also, compute the initial guess.
99: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
101: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
102: Create nonlinear solver context
104: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
105: DMSetApplicationContext(da,&user);
106: DMDASNESSetFunctionLocal(da,INSERT_VALUES,(PetscErrorCode (*)(DMDALocalInfo*,void*,void*,void*))FormFunctionLocal,&user);
108: if (errorinmatmult) {
109: MatCreateSNESMF(user.snes,&Jmf);
110: MatSetFromOptions(Jmf);
111: MatGetLocalSize(Jmf,&mlocal,&nlocal);
112: matshellctx.Jmf = Jmf;
113: MatCreateShell(PetscObjectComm((PetscObject)Jmf),mlocal,nlocal,PETSC_DECIDE,PETSC_DECIDE,&matshellctx,&J);
114: MatShellSetOperation(J,MATOP_MULT,(void (*)(void))MatMult_MyShell);
115: MatShellSetOperation(J,MATOP_ASSEMBLY_END,(void (*)(void))MatAssemblyEnd_MyShell);
116: SNESSetJacobian(user.snes,J,J,MatMFFDComputeJacobian,NULL);
117: }
119: SNESSetFromOptions(user.snes);
120: PetscPrintf(comm,"lid velocity = %g, prandtl # = %g, grashof # = %g\n",(double)user.lidvelocity,(double)user.prandtl,(double)user.grashof);
122: if (errorinpcapply) {
123: SNESGetKSP(user.snes,&ksp);
124: KSPGetPC(ksp,&pc);
125: PCSetType(pc,PCSHELL);
126: PCShellSetApply(pc,PCApply_MyShell);
127: }
129: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
130: Solve the nonlinear system
131: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
132: DMCreateGlobalVector(da,&x);
133: FormInitialGuess(&user,da,x);
135: if (errorinpcsetup) {
136: SNESSetUp(user.snes);
137: SNESSetJacobian(user.snes,NULL,NULL,SNESComputeJacobian_MyShell,NULL);
138: }
139: SNESSolve(user.snes,NULL,x);
142: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
143: Free work space. All PETSc objects should be destroyed when they
144: are no longer needed.
145: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
146: MatDestroy(&J);
147: MatDestroy(&Jmf);
148: VecDestroy(&x);
149: DMDestroy(&da);
150: SNESDestroy(&user.snes);
151: PetscFinalize();
152: return 0;
153: }
155: /* ------------------------------------------------------------------- */
160: /*
161: FormInitialGuess - Forms initial approximation.
163: Input Parameters:
164: user - user-defined application context
165: X - vector
167: Output Parameter:
168: X - vector
169: */
170: PetscErrorCode FormInitialGuess(AppCtx *user,DM da,Vec X)
171: {
172: PetscInt i,j,mx,xs,ys,xm,ym;
174: PetscReal grashof,dx;
175: Field **x;
178: grashof = user->grashof;
180: DMDAGetInfo(da,0,&mx,0,0,0,0,0,0,0,0,0,0,0);
181: dx = 1.0/(mx-1);
183: /*
184: Get local grid boundaries (for 2-dimensional DMDA):
185: xs, ys - starting grid indices (no ghost points)
186: xm, ym - widths of local grid (no ghost points)
187: */
188: DMDAGetCorners(da,&xs,&ys,NULL,&xm,&ym,NULL);
190: /*
191: Get a pointer to vector data.
192: - For default PETSc vectors, VecGetArray() returns a pointer to
193: the data array. Otherwise, the routine is implementation dependent.
194: - You MUST call VecRestoreArray() when you no longer need access to
195: the array.
196: */
197: DMDAVecGetArray(da,X,&x);
199: /*
200: Compute initial guess over the locally owned part of the grid
201: Initial condition is motionless fluid and equilibrium temperature
202: */
203: for (j=ys; j<ys+ym; j++) {
204: for (i=xs; i<xs+xm; i++) {
205: x[j][i].u = 0.0;
206: x[j][i].v = 0.0;
207: x[j][i].omega = 0.0;
208: x[j][i].temp = (grashof>0)*i*dx;
209: }
210: }
212: /*
213: Restore vector
214: */
215: DMDAVecRestoreArray(da,X,&x);
216: return(0);
217: }
221: PetscErrorCode FormFunctionLocal(DMDALocalInfo *info,Field **x,Field **f,void *ptr)
222: {
223: AppCtx *user = (AppCtx*)ptr;
224: PetscErrorCode ierr;
225: PetscInt xints,xinte,yints,yinte,i,j;
226: PetscReal hx,hy,dhx,dhy,hxdhy,hydhx;
227: PetscReal grashof,prandtl,lid;
228: PetscScalar u,uxx,uyy,vx,vy,avx,avy,vxp,vxm,vyp,vym;
229: static PetscInt fail = 0;
232: if ((fail++ > 7 && user->errorindomainmf) || (fail++ > 36 && user->errorindomain)){
233: PetscMPIInt rank;
234: MPI_Comm_rank(PetscObjectComm((PetscObject)user->snes),&rank);
235: if (!rank) {
236: SNESSetFunctionDomainError(user->snes);
237: }
238: }
239: grashof = user->grashof;
240: prandtl = user->prandtl;
241: lid = user->lidvelocity;
243: /*
244: Define mesh intervals ratios for uniform grid.
246: Note: FD formulae below are normalized by multiplying through by
247: local volume element (i.e. hx*hy) to obtain coefficients O(1) in two dimensions.
250: */
251: dhx = (PetscReal)(info->mx-1); dhy = (PetscReal)(info->my-1);
252: hx = 1.0/dhx; hy = 1.0/dhy;
253: hxdhy = hx*dhy; hydhx = hy*dhx;
255: xints = info->xs; xinte = info->xs+info->xm; yints = info->ys; yinte = info->ys+info->ym;
257: /* Test whether we are on the bottom edge of the global array */
258: if (yints == 0) {
259: j = 0;
260: yints = yints + 1;
261: /* bottom edge */
262: for (i=info->xs; i<info->xs+info->xm; i++) {
263: f[j][i].u = x[j][i].u;
264: f[j][i].v = x[j][i].v;
265: f[j][i].omega = x[j][i].omega + (x[j+1][i].u - x[j][i].u)*dhy;
266: f[j][i].temp = x[j][i].temp-x[j+1][i].temp;
267: }
268: }
270: /* Test whether we are on the top edge of the global array */
271: if (yinte == info->my) {
272: j = info->my - 1;
273: yinte = yinte - 1;
274: /* top edge */
275: for (i=info->xs; i<info->xs+info->xm; i++) {
276: f[j][i].u = x[j][i].u - lid;
277: f[j][i].v = x[j][i].v;
278: f[j][i].omega = x[j][i].omega + (x[j][i].u - x[j-1][i].u)*dhy;
279: f[j][i].temp = x[j][i].temp-x[j-1][i].temp;
280: }
281: }
283: /* Test whether we are on the left edge of the global array */
284: if (xints == 0) {
285: i = 0;
286: xints = xints + 1;
287: /* left edge */
288: for (j=info->ys; j<info->ys+info->ym; j++) {
289: f[j][i].u = x[j][i].u;
290: f[j][i].v = x[j][i].v;
291: f[j][i].omega = x[j][i].omega - (x[j][i+1].v - x[j][i].v)*dhx;
292: f[j][i].temp = x[j][i].temp;
293: }
294: }
296: /* Test whether we are on the right edge of the global array */
297: if (xinte == info->mx) {
298: i = info->mx - 1;
299: xinte = xinte - 1;
300: /* right edge */
301: for (j=info->ys; j<info->ys+info->ym; j++) {
302: f[j][i].u = x[j][i].u;
303: f[j][i].v = x[j][i].v;
304: f[j][i].omega = x[j][i].omega - (x[j][i].v - x[j][i-1].v)*dhx;
305: f[j][i].temp = x[j][i].temp - (PetscReal)(grashof>0);
306: }
307: }
309: /* Compute over the interior points */
310: for (j=yints; j<yinte; j++) {
311: for (i=xints; i<xinte; i++) {
313: /*
314: convective coefficients for upwinding
315: */
316: vx = x[j][i].u; avx = PetscAbsScalar(vx);
317: vxp = .5*(vx+avx); vxm = .5*(vx-avx);
318: vy = x[j][i].v; avy = PetscAbsScalar(vy);
319: vyp = .5*(vy+avy); vym = .5*(vy-avy);
321: /* U velocity */
322: u = x[j][i].u;
323: uxx = (2.0*u - x[j][i-1].u - x[j][i+1].u)*hydhx;
324: uyy = (2.0*u - x[j-1][i].u - x[j+1][i].u)*hxdhy;
325: f[j][i].u = uxx + uyy - .5*(x[j+1][i].omega-x[j-1][i].omega)*hx;
327: /* V velocity */
328: u = x[j][i].v;
329: uxx = (2.0*u - x[j][i-1].v - x[j][i+1].v)*hydhx;
330: uyy = (2.0*u - x[j-1][i].v - x[j+1][i].v)*hxdhy;
331: f[j][i].v = uxx + uyy + .5*(x[j][i+1].omega-x[j][i-1].omega)*hy;
333: /* Omega */
334: u = x[j][i].omega;
335: uxx = (2.0*u - x[j][i-1].omega - x[j][i+1].omega)*hydhx;
336: uyy = (2.0*u - x[j-1][i].omega - x[j+1][i].omega)*hxdhy;
337: f[j][i].omega = uxx + uyy + (vxp*(u - x[j][i-1].omega) + vxm*(x[j][i+1].omega - u))*hy +
338: (vyp*(u - x[j-1][i].omega) + vym*(x[j+1][i].omega - u))*hx -
339: .5*grashof*(x[j][i+1].temp - x[j][i-1].temp)*hy;
341: /* Temperature */
342: u = x[j][i].temp;
343: uxx = (2.0*u - x[j][i-1].temp - x[j][i+1].temp)*hydhx;
344: uyy = (2.0*u - x[j-1][i].temp - x[j+1][i].temp)*hxdhy;
345: f[j][i].temp = uxx + uyy + prandtl*((vxp*(u - x[j][i-1].temp) + vxm*(x[j][i+1].temp - u))*hy +
346: (vyp*(u - x[j-1][i].temp) + vym*(x[j+1][i].temp - u))*hx);
347: }
348: }
350: /*
351: Flop count (multiply-adds are counted as 2 operations)
352: */
353: PetscLogFlops(84.0*info->ym*info->xm);
354: return(0);
355: }
359: PetscErrorCode MatMult_MyShell(Mat A,Vec x,Vec y)
360: {
361: PetscErrorCode ierr;
362: MatShellCtx *matshellctx;
363: static PetscInt fail = 0;
366: MatShellGetContext(A,&matshellctx);
367: MatMult(matshellctx->Jmf,x,y);
368: if (fail++ > 5) {
369: PetscMPIInt rank;
370: MPI_Comm_rank(PetscObjectComm((PetscObject)A),&rank);
371: if (!rank) {VecSetInf(y);}
372: }
373: return(0);
374: }
378: PetscErrorCode MatAssemblyEnd_MyShell(Mat A,MatAssemblyType tp)
379: {
381: MatShellCtx *matshellctx;
384: MatShellGetContext(A,&matshellctx);
385: MatAssemblyEnd(matshellctx->Jmf,tp);
386: return(0);
387: }
391: PetscErrorCode PCApply_MyShell(PC pc,Vec x,Vec y)
392: {
394: static PetscInt fail = 0;
397: VecCopy(x,y);
398: if (fail++ > 3) {
399: PetscMPIInt rank;
400: MPI_Comm_rank(PetscObjectComm((PetscObject)pc),&rank);
401: if (!rank) {VecSetInf(y);}
402: }
403: return(0);
404: }
406: PETSC_EXTERN PetscErrorCode SNESComputeJacobian_DMDA(SNES,Vec,Mat,Mat,void*);
410: PetscErrorCode SNESComputeJacobian_MyShell(SNES snes,Vec X,Mat A,Mat B,void *ctx)
411: {
412: static PetscInt fail = 0;
413: PetscErrorCode ierr;
416: SNESComputeJacobian_DMDA(snes,X,A,B,ctx);
417: if (fail++ > 0) {
418: MatZeroEntries(A);
419: }
420: return(0);
421: }