12    User-defined routines and data structures
22   PetscBool draw_contours;                 /* flag - 1 indicates drawing contours */
40   AppCtx      user; /* user-defined work context */  in main()
45   Mat         J = NULL, Jmf = NULL;  in main()  local
54   PetscCall(PetscOptionsGetBool(NULL, NULL, "-error_in_matmult", &errorinmatmult, NULL));  in main()
55   PetscCall(PetscOptionsGetBool(NULL, NULL, "-error_in_pcapply", &errorinpcapply, NULL));  in main()
56   PetscCall(PetscOptionsGetBool(NULL, NULL, "-error_in_pcsetup", &errorinpcsetup, NULL));  in main()
58   PetscCall(PetscOptionsGetBool(NULL, NULL, "-error_in_domain", &user.errorindomain, NULL));  in main()
60   PetscCall(PetscOptionsGetBool(NULL, NULL, "-error_in_domainmf", &user.errorindomainmf, NULL));  in main()
82   PetscCall(PetscOptionsGetReal(NULL, NULL, "-lidvelocity", &user.lidvelocity, NULL));  in main()
83   PetscCall(PetscOptionsGetReal(NULL, NULL, "-prandtl", &user.prandtl, NULL));  in main()
84   PetscCall(PetscOptionsGetReal(NULL, NULL, "-grashof", &user.grashof, NULL));  in main()
85   PetscCall(PetscOptionsHasName(NULL, NULL, "-contours", &user.draw_contours));  in main()
92   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -  in main()
95      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */  in main()
97   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -  in main()
100      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */  in main()
109 …(PetscObjectComm((PetscObject)Jmf), mlocal, nlocal, PETSC_DECIDE, PETSC_DECIDE, &matshellctx, &J));  in main()
110     PetscCall(MatShellSetOperation(J, MATOP_MULT, (PetscErrorCodeFn *)MatMult_MyShell));  in main()
111 …PetscCall(MatShellSetOperation(J, MATOP_ASSEMBLY_END, (PetscErrorCodeFn *)MatAssemblyEnd_MyShell));  in main()
112     PetscCall(SNESSetJacobian(user.snes, J, J, MatMFFDComputeJacobian, NULL));  in main()
125   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -  in main()
127      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */  in main()
137   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -  in main()
140      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */  in main()
141   PetscCall(MatDestroy(&J));  in main()
151    FormInitialGuess - Forms initial approximation.
154    user - user-defined application context
155    X - vector
158    X - vector
162   PetscInt  i, j, mx, xs, ys, xm, ym;  in FormInitialGuess()  local
167   grashof = user->grashof;  in FormInitialGuess()
170   dx = 1.0 / (mx - 1);  in FormInitialGuess()
173      Get local grid boundaries (for 2-dimensional DMDA):  in FormInitialGuess()
174        xs, ys   - starting grid indices (no ghost points)  in FormInitialGuess()
175        xm, ym   - widths of local grid (no ghost points)  in FormInitialGuess()
181        - For default PETSc vectors, VecGetArray() returns a pointer to  in FormInitialGuess()
183        - You MUST call VecRestoreArray() when you no longer need access to  in FormInitialGuess()
192   for (j = ys; j < ys + ym; j++) {  in FormInitialGuess()
194       x[j][i].u     = 0.0;  in FormInitialGuess()
195       x[j][i].v     = 0.0;  in FormInitialGuess()
196       x[j][i].omega = 0.0;  in FormInitialGuess()
197       x[j][i].temp  = (grashof > 0) * i * dx;  in FormInitialGuess()
211   PetscInt        xints, xinte, yints, yinte, i, j;  in FormFunctionLocal()  local
218   if ((fail++ > 7 && user->errorindomainmf) || (fail++ > 36 && user->errorindomain)) {  in FormFunctionLocal()
220     PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)user->snes), &rank));  in FormFunctionLocal()
221     if (rank == 0) PetscCall(SNESSetFunctionDomainError(user->snes));  in FormFunctionLocal()
223   grashof = user->grashof;  in FormFunctionLocal()
224   prandtl = user->prandtl;  in FormFunctionLocal()
225   lid     = user->lidvelocity;  in FormFunctionLocal()
234   dhx   = (PetscReal)(info->mx - 1);  in FormFunctionLocal()
235   dhy   = (PetscReal)(info->my - 1);  in FormFunctionLocal()
241   xints = info->xs;  in FormFunctionLocal()
242   xinte = info->xs + info->xm;  in FormFunctionLocal()
243   yints = info->ys;  in FormFunctionLocal()
244   yinte = info->ys + info->ym;  in FormFunctionLocal()
248     j     = 0;  in FormFunctionLocal()
251     for (i = info->xs; i < info->xs + info->xm; i++) {  in FormFunctionLocal()
252       f[j][i].u     = x[j][i].u;  in FormFunctionLocal()
253       f[j][i].v     = x[j][i].v;  in FormFunctionLocal()
254       f[j][i].omega = x[j][i].omega + (x[j + 1][i].u - x[j][i].u) * dhy;  in FormFunctionLocal()
255       f[j][i].temp  = x[j][i].temp - x[j + 1][i].temp;  in FormFunctionLocal()
260   if (yinte == info->my) {  in FormFunctionLocal()
261     j     = info->my - 1;  in FormFunctionLocal()
262     yinte = yinte - 1;  in FormFunctionLocal()
264     for (i = info->xs; i < info->xs + info->xm; i++) {  in FormFunctionLocal()
265       f[j][i].u     = x[j][i].u - lid;  in FormFunctionLocal()
266       f[j][i].v     = x[j][i].v;  in FormFunctionLocal()
267       f[j][i].omega = x[j][i].omega + (x[j][i].u - x[j - 1][i].u) * dhy;  in FormFunctionLocal()
268       f[j][i].temp  = x[j][i].temp - x[j - 1][i].temp;  in FormFunctionLocal()
277     for (j = info->ys; j < info->ys + info->ym; j++) {  in FormFunctionLocal()
278       f[j][i].u     = x[j][i].u;  in FormFunctionLocal()
279       f[j][i].v     = x[j][i].v;  in FormFunctionLocal()
280       f[j][i].omega = x[j][i].omega - (x[j][i + 1].v - x[j][i].v) * dhx;  in FormFunctionLocal()
281       f[j][i].temp  = x[j][i].temp;  in FormFunctionLocal()
286   if (xinte == info->mx) {  in FormFunctionLocal()
287     i     = info->mx - 1;  in FormFunctionLocal()
288     xinte = xinte - 1;  in FormFunctionLocal()
290     for (j = info->ys; j < info->ys + info->ym; j++) {  in FormFunctionLocal()
291       f[j][i].u     = x[j][i].u;  in FormFunctionLocal()
292       f[j][i].v     = x[j][i].v;  in FormFunctionLocal()
293       f[j][i].omega = x[j][i].omega - (x[j][i].v - x[j][i - 1].v) * dhx;  in FormFunctionLocal()
294       f[j][i].temp  = x[j][i].temp - (PetscReal)(grashof > 0);  in FormFunctionLocal()
299   for (j = yints; j < yinte; j++) {  in FormFunctionLocal()
304       vx  = x[j][i].u;  in FormFunctionLocal()
307       vxm = .5 * (vx - avx);  in FormFunctionLocal()
308       vy  = x[j][i].v;  in FormFunctionLocal()
311       vym = .5 * (vy - avy);  in FormFunctionLocal()
314       u         = x[j][i].u;  in FormFunctionLocal()
315       uxx       = (2.0 * u - x[j][i - 1].u - x[j][i + 1].u) * hydhx;  in FormFunctionLocal()
316       uyy       = (2.0 * u - x[j - 1][i].u - x[j + 1][i].u) * hxdhy;  in FormFunctionLocal()
317       f[j][i].u = uxx + uyy - .5 * (x[j + 1][i].omega - x[j - 1][i].omega) * hx;  in FormFunctionLocal()
320       u         = x[j][i].v;  in FormFunctionLocal()
321       uxx       = (2.0 * u - x[j][i - 1].v - x[j][i + 1].v) * hydhx;  in FormFunctionLocal()
322       uyy       = (2.0 * u - x[j - 1][i].v - x[j + 1][i].v) * hxdhy;  in FormFunctionLocal()
323       f[j][i].v = uxx + uyy + .5 * (x[j][i + 1].omega - x[j][i - 1].omega) * hy;  in FormFunctionLocal()
326       u             = x[j][i].omega;  in FormFunctionLocal()
327       uxx           = (2.0 * u - x[j][i - 1].omega - x[j][i + 1].omega) * hydhx;  in FormFunctionLocal()
328       uyy           = (2.0 * u - x[j - 1][i].omega - x[j + 1][i].omega) * hxdhy;  in FormFunctionLocal()
329 …j][i].omega = uxx + uyy + (vxp * (u - x[j][i - 1].omega) + vxm * (x[j][i + 1].omega - u)) * hy + (…  in FormFunctionLocal()
332       u            = x[j][i].temp;  in FormFunctionLocal()
333       uxx          = (2.0 * u - x[j][i - 1].temp - x[j][i + 1].temp) * hydhx;  in FormFunctionLocal()
334       uyy          = (2.0 * u - x[j - 1][i].temp - x[j + 1][i].temp) * hxdhy;  in FormFunctionLocal()
335 …j][i].temp = uxx + uyy + prandtl * ((vxp * (u - x[j][i - 1].temp) + vxm * (x[j][i + 1].temp - u)) …  in FormFunctionLocal()
340      Flop count (multiply-adds are counted as 2 operations)  in FormFunctionLocal()
342   PetscCall(PetscLogFlops(84.0 * info->ym * info->xm));  in FormFunctionLocal()
354   PetscCall(MatMult(matshellctx->Jmf, x, y));  in MatMult_MyShell()
370   PetscCall(MatAssemblyEnd(matshellctx->Jmf, tp));  in MatAssemblyEnd_MyShell()
405       args: -snes_converged_reason -ksp_converged_reason
409       args: -snes_converged_reason -ksp_converged_reason -error_in_matmult -fp_trap 0
413       args: -snes_converged_reason -ksp_converged_reason -error_in_pcapply -fp_trap 0
417       args: -snes_converged_reason -ksp_converged_reason -error_in_pcsetup -fp_trap 0
421 …   args: -snes_converged_reason -ksp_converged_reason -error_in_pcsetup -pc_type bjacobi -fp_trap 0
425 …args: -snes_converged_reason -ksp_converged_reason -error_in_pcsetup -pc_type fieldsplit -fp_trap 0
430 …args: -snes_converged_reason -ksp_converged_reason -error_in_domainmf -snes_mf -pc_type none -fp_t…
434       args: -snes_converged_reason -ksp_converged_reason -error_in_domain -fp_trap 0
438       args: -snes_converged_reason -ksp_converged_reason -error_in_domain -snes_mf -pc_type none