Actual source code: ex69.c
1: static char help[] = "Tests recovery from domain errors in MatMult() and PCApply()\n\n";
3: /*
4: See src/ksp/ksp/tutorials/ex19.c from which this was copied
5: */
7: #include <petscsnes.h>
8: #include <petscdm.h>
9: #include <petscdmda.h>
11: /*
12: User-defined routines and data structures
13: */
14: typedef struct {
15: PetscScalar u, v, omega, temp;
16: } Field;
18: PetscErrorCode FormFunctionLocal(DMDALocalInfo *, Field **, Field **, void *);
20: typedef struct {
21: PetscReal lidvelocity, prandtl, grashof; /* physical parameters */
22: PetscBool draw_contours; /* flag - 1 indicates drawing contours */
23: PetscBool errorindomain;
24: PetscBool errorindomainmf;
25: SNES snes;
26: } AppCtx;
28: typedef struct {
29: Mat Jmf;
30: } MatShellCtx;
32: extern PetscErrorCode FormInitialGuess(AppCtx *, DM, Vec);
33: extern PetscErrorCode MatMult_MyShell(Mat, Vec, Vec);
34: extern PetscErrorCode MatAssemblyEnd_MyShell(Mat, MatAssemblyType);
35: extern PetscErrorCode PCApply_MyShell(PC, Vec, Vec);
36: extern PetscErrorCode SNESComputeJacobian_MyShell(SNES, Vec, Mat, Mat, void *);
38: int main(int argc, char **argv)
39: {
40: AppCtx user; /* user-defined work context */
41: PetscInt mx, my;
42: MPI_Comm comm;
43: DM da;
44: Vec x;
45: Mat J = NULL, Jmf = NULL;
46: MatShellCtx matshellctx;
47: PetscInt mlocal, nlocal;
48: PC pc;
49: KSP ksp;
50: PetscBool errorinmatmult = PETSC_FALSE, errorinpcapply = PETSC_FALSE, errorinpcsetup = PETSC_FALSE;
52: PetscFunctionBeginUser;
53: PetscCall(PetscInitialize(&argc, &argv, NULL, help));
54: PetscCall(PetscOptionsGetBool(NULL, NULL, "-error_in_matmult", &errorinmatmult, NULL));
55: PetscCall(PetscOptionsGetBool(NULL, NULL, "-error_in_pcapply", &errorinpcapply, NULL));
56: PetscCall(PetscOptionsGetBool(NULL, NULL, "-error_in_pcsetup", &errorinpcsetup, NULL));
57: user.errorindomain = PETSC_FALSE;
58: PetscCall(PetscOptionsGetBool(NULL, NULL, "-error_in_domain", &user.errorindomain, NULL));
59: user.errorindomainmf = PETSC_FALSE;
60: PetscCall(PetscOptionsGetBool(NULL, NULL, "-error_in_domainmf", &user.errorindomainmf, NULL));
62: comm = PETSC_COMM_WORLD;
63: PetscCall(SNESCreate(comm, &user.snes));
65: /*
66: Create distributed array object to manage parallel grid and vectors
67: for principal unknowns (x) and governing residuals (f)
68: */
69: PetscCall(DMDACreate2d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DMDA_STENCIL_STAR, 4, 4, PETSC_DECIDE, PETSC_DECIDE, 4, 1, 0, 0, &da));
70: PetscCall(DMSetFromOptions(da));
71: PetscCall(DMSetUp(da));
72: PetscCall(SNESSetDM(user.snes, da));
74: PetscCall(DMDAGetInfo(da, 0, &mx, &my, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE));
75: /*
76: Problem parameters (velocity of lid, prandtl, and grashof numbers)
77: */
78: user.lidvelocity = 1.0 / (mx * my);
79: user.prandtl = 1.0;
80: user.grashof = 1.0;
82: PetscCall(PetscOptionsGetReal(NULL, NULL, "-lidvelocity", &user.lidvelocity, NULL));
83: PetscCall(PetscOptionsGetReal(NULL, NULL, "-prandtl", &user.prandtl, NULL));
84: PetscCall(PetscOptionsGetReal(NULL, NULL, "-grashof", &user.grashof, NULL));
85: PetscCall(PetscOptionsHasName(NULL, NULL, "-contours", &user.draw_contours));
87: PetscCall(DMDASetFieldName(da, 0, "x_velocity"));
88: PetscCall(DMDASetFieldName(da, 1, "y_velocity"));
89: PetscCall(DMDASetFieldName(da, 2, "Omega"));
90: PetscCall(DMDASetFieldName(da, 3, "temperature"));
92: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
93: Create user context, set problem data, create vector data structures.
94: Also, compute the initial guess.
95: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
97: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
98: Create nonlinear solver context
100: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
101: PetscCall(DMSetApplicationContext(da, &user));
102: PetscCall(DMDASNESSetFunctionLocal(da, INSERT_VALUES, (PetscErrorCode (*)(DMDALocalInfo *, void *, void *, void *))FormFunctionLocal, &user));
104: if (errorinmatmult) {
105: PetscCall(MatCreateSNESMF(user.snes, &Jmf));
106: PetscCall(MatSetFromOptions(Jmf));
107: PetscCall(MatGetLocalSize(Jmf, &mlocal, &nlocal));
108: matshellctx.Jmf = Jmf;
109: PetscCall(MatCreateShell(PetscObjectComm((PetscObject)Jmf), mlocal, nlocal, PETSC_DECIDE, PETSC_DECIDE, &matshellctx, &J));
110: PetscCall(MatShellSetOperation(J, MATOP_MULT, (void (*)(void))MatMult_MyShell));
111: PetscCall(MatShellSetOperation(J, MATOP_ASSEMBLY_END, (void (*)(void))MatAssemblyEnd_MyShell));
112: PetscCall(SNESSetJacobian(user.snes, J, J, MatMFFDComputeJacobian, NULL));
113: }
115: PetscCall(SNESSetFromOptions(user.snes));
116: PetscCall(PetscPrintf(comm, "lid velocity = %g, prandtl # = %g, grashof # = %g\n", (double)user.lidvelocity, (double)user.prandtl, (double)user.grashof));
118: if (errorinpcapply) {
119: PetscCall(SNESGetKSP(user.snes, &ksp));
120: PetscCall(KSPGetPC(ksp, &pc));
121: PetscCall(PCSetType(pc, PCSHELL));
122: PetscCall(PCShellSetApply(pc, PCApply_MyShell));
123: }
125: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
126: Solve the nonlinear system
127: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
128: PetscCall(DMCreateGlobalVector(da, &x));
129: PetscCall(FormInitialGuess(&user, da, x));
131: if (errorinpcsetup) {
132: PetscCall(SNESSetUp(user.snes));
133: PetscCall(SNESSetJacobian(user.snes, NULL, NULL, SNESComputeJacobian_MyShell, NULL));
134: }
135: PetscCall(SNESSolve(user.snes, NULL, x));
137: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
138: Free work space. All PETSc objects should be destroyed when they
139: are no longer needed.
140: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
141: PetscCall(MatDestroy(&J));
142: PetscCall(MatDestroy(&Jmf));
143: PetscCall(VecDestroy(&x));
144: PetscCall(DMDestroy(&da));
145: PetscCall(SNESDestroy(&user.snes));
146: PetscCall(PetscFinalize());
147: return 0;
148: }
150: /*
151: FormInitialGuess - Forms initial approximation.
153: Input Parameters:
154: user - user-defined application context
155: X - vector
157: Output Parameter:
158: X - vector
159: */
160: PetscErrorCode FormInitialGuess(AppCtx *user, DM da, Vec X)
161: {
162: PetscInt i, j, mx, xs, ys, xm, ym;
163: PetscReal grashof, dx;
164: Field **x;
166: PetscFunctionBeginUser;
167: grashof = user->grashof;
169: PetscCall(DMDAGetInfo(da, 0, &mx, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0));
170: dx = 1.0 / (mx - 1);
172: /*
173: Get local grid boundaries (for 2-dimensional DMDA):
174: xs, ys - starting grid indices (no ghost points)
175: xm, ym - widths of local grid (no ghost points)
176: */
177: PetscCall(DMDAGetCorners(da, &xs, &ys, NULL, &xm, &ym, NULL));
179: /*
180: Get a pointer to vector data.
181: - For default PETSc vectors, VecGetArray() returns a pointer to
182: the data array. Otherwise, the routine is implementation dependent.
183: - You MUST call VecRestoreArray() when you no longer need access to
184: the array.
185: */
186: PetscCall(DMDAVecGetArray(da, X, &x));
188: /*
189: Compute initial guess over the locally owned part of the grid
190: Initial condition is motionless fluid and equilibrium temperature
191: */
192: for (j = ys; j < ys + ym; j++) {
193: for (i = xs; i < xs + xm; i++) {
194: x[j][i].u = 0.0;
195: x[j][i].v = 0.0;
196: x[j][i].omega = 0.0;
197: x[j][i].temp = (grashof > 0) * i * dx;
198: }
199: }
201: /*
202: Restore vector
203: */
204: PetscCall(DMDAVecRestoreArray(da, X, &x));
205: PetscFunctionReturn(PETSC_SUCCESS);
206: }
208: PetscErrorCode FormFunctionLocal(DMDALocalInfo *info, Field **x, Field **f, void *ptr)
209: {
210: AppCtx *user = (AppCtx *)ptr;
211: PetscInt xints, xinte, yints, yinte, i, j;
212: PetscReal hx, hy, dhx, dhy, hxdhy, hydhx;
213: PetscReal grashof, prandtl, lid;
214: PetscScalar u, uxx, uyy, vx, vy, avx, avy, vxp, vxm, vyp, vym;
215: static PetscInt fail = 0;
217: PetscFunctionBeginUser;
218: if ((fail++ > 7 && user->errorindomainmf) || (fail++ > 36 && user->errorindomain)) {
219: PetscMPIInt rank;
220: PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)user->snes), &rank));
221: if (rank == 0) PetscCall(SNESSetFunctionDomainError(user->snes));
222: }
223: grashof = user->grashof;
224: prandtl = user->prandtl;
225: lid = user->lidvelocity;
227: /*
228: Define mesh intervals ratios for uniform grid.
230: Note: FD formulae below are normalized by multiplying through by
231: local volume element (i.e. hx*hy) to obtain coefficients O(1) in two dimensions.
233: */
234: dhx = (PetscReal)(info->mx - 1);
235: dhy = (PetscReal)(info->my - 1);
236: hx = 1.0 / dhx;
237: hy = 1.0 / dhy;
238: hxdhy = hx * dhy;
239: hydhx = hy * dhx;
241: xints = info->xs;
242: xinte = info->xs + info->xm;
243: yints = info->ys;
244: 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++) {
301: /*
302: convective coefficients for upwinding
303: */
304: vx = x[j][i].u;
305: avx = PetscAbsScalar(vx);
306: vxp = .5 * (vx + avx);
307: vxm = .5 * (vx - avx);
308: vy = x[j][i].v;
309: avy = PetscAbsScalar(vy);
310: vyp = .5 * (vy + avy);
311: vym = .5 * (vy - avy);
313: /* U velocity */
314: u = x[j][i].u;
315: uxx = (2.0 * u - x[j][i - 1].u - x[j][i + 1].u) * hydhx;
316: uyy = (2.0 * u - x[j - 1][i].u - x[j + 1][i].u) * hxdhy;
317: f[j][i].u = uxx + uyy - .5 * (x[j + 1][i].omega - x[j - 1][i].omega) * hx;
319: /* V velocity */
320: u = x[j][i].v;
321: uxx = (2.0 * u - x[j][i - 1].v - x[j][i + 1].v) * hydhx;
322: uyy = (2.0 * u - x[j - 1][i].v - x[j + 1][i].v) * hxdhy;
323: f[j][i].v = uxx + uyy + .5 * (x[j][i + 1].omega - x[j][i - 1].omega) * hy;
325: /* Omega */
326: u = x[j][i].omega;
327: uxx = (2.0 * u - x[j][i - 1].omega - x[j][i + 1].omega) * hydhx;
328: uyy = (2.0 * u - x[j - 1][i].omega - x[j + 1][i].omega) * hxdhy;
329: f[j][i].omega = uxx + uyy + (vxp * (u - x[j][i - 1].omega) + vxm * (x[j][i + 1].omega - u)) * hy + (vyp * (u - x[j - 1][i].omega) + vym * (x[j + 1][i].omega - u)) * hx - .5 * grashof * (x[j][i + 1].temp - x[j][i - 1].temp) * hy;
331: /* Temperature */
332: u = x[j][i].temp;
333: uxx = (2.0 * u - x[j][i - 1].temp - x[j][i + 1].temp) * hydhx;
334: uyy = (2.0 * u - x[j - 1][i].temp - x[j + 1][i].temp) * hxdhy;
335: f[j][i].temp = uxx + uyy + prandtl * ((vxp * (u - x[j][i - 1].temp) + vxm * (x[j][i + 1].temp - u)) * hy + (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: PetscCall(PetscLogFlops(84.0 * info->ym * info->xm));
343: PetscFunctionReturn(PETSC_SUCCESS);
344: }
346: PetscErrorCode MatMult_MyShell(Mat A, Vec x, Vec y)
347: {
348: MatShellCtx *matshellctx;
349: static PetscInt fail = 0;
350: PetscMPIInt rank;
352: PetscFunctionBegin;
353: PetscCall(MatShellGetContext(A, &matshellctx));
354: PetscCall(MatMult(matshellctx->Jmf, x, y));
355: PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)A), &rank));
356: if (fail++ > 5) {
357: PetscCall(VecFlag(y, rank == 0));
358: PetscCall(VecAssemblyBegin(y));
359: PetscCall(VecAssemblyEnd(y));
360: }
361: PetscFunctionReturn(PETSC_SUCCESS);
362: }
364: PetscErrorCode MatAssemblyEnd_MyShell(Mat A, MatAssemblyType tp)
365: {
366: MatShellCtx *matshellctx;
368: PetscFunctionBegin;
369: PetscCall(MatShellGetContext(A, &matshellctx));
370: PetscCall(MatAssemblyEnd(matshellctx->Jmf, tp));
371: PetscFunctionReturn(PETSC_SUCCESS);
372: }
374: PetscErrorCode PCApply_MyShell(PC pc, Vec x, Vec y)
375: {
376: static PetscInt fail = 0;
377: PetscMPIInt rank;
379: PetscFunctionBegin;
380: PetscCall(VecCopy(x, y));
381: PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)pc), &rank));
382: if (fail++ > 3) {
383: PetscCall(VecFlag(y, rank == 0));
384: PetscCall(VecAssemblyBegin(y));
385: PetscCall(VecAssemblyEnd(y));
386: }
387: PetscFunctionReturn(PETSC_SUCCESS);
388: }
390: PETSC_EXTERN PetscErrorCode SNESComputeJacobian_DMDA(SNES, Vec, Mat, Mat, void *);
392: PetscErrorCode SNESComputeJacobian_MyShell(SNES snes, Vec X, Mat A, Mat B, void *ctx)
393: {
394: static PetscInt fail = 0;
396: PetscFunctionBegin;
397: PetscCall(SNESComputeJacobian_DMDA(snes, X, A, B, ctx));
398: if (fail++ > 0) PetscCall(MatZeroEntries(A));
399: PetscFunctionReturn(PETSC_SUCCESS);
400: }
402: /*TEST
404: test:
405: args: -snes_converged_reason -ksp_converged_reason
407: test:
408: suffix: 2
409: args: -snes_converged_reason -ksp_converged_reason -error_in_matmult -fp_trap 0
411: test:
412: suffix: 3
413: args: -snes_converged_reason -ksp_converged_reason -error_in_pcapply -fp_trap 0
415: test:
416: suffix: 4
417: args: -snes_converged_reason -ksp_converged_reason -error_in_pcsetup -fp_trap 0
419: test:
420: suffix: 5
421: args: -snes_converged_reason -ksp_converged_reason -error_in_pcsetup -pc_type bjacobi -fp_trap 0
423: test:
424: suffix: 5_fieldsplit
425: args: -snes_converged_reason -ksp_converged_reason -error_in_pcsetup -pc_type fieldsplit -fp_trap 0
426: output_file: output/ex69_5.out
428: test:
429: suffix: 6
430: args: -snes_converged_reason -ksp_converged_reason -error_in_domainmf -snes_mf -pc_type none -fp_trap 0
432: test:
433: suffix: 7
434: args: -snes_converged_reason -ksp_converged_reason -error_in_domain -fp_trap 0
436: test:
437: suffix: 8
438: args: -snes_converged_reason -ksp_converged_reason -error_in_domain -snes_mf -pc_type none
440: TEST*/