v0.14.0
Loading...
Searching...
No Matches
HackMG.cpp
Go to the documentation of this file.
1#include <MoFEM.hpp>
2
3#if PETSC_VERSION_GE(3, 6, 0)
4#include <petsc/private/kspimpl.h> /*I "petscksp.h" I*/
5#include <petsc/private/dmimpl.h> /*I "petscdm.h" I*/
6#include <petsc/private/pcmgimpl.h>
7#else
8#error "PETSc version should be 3.6.0 or higher"
9#endif
10
12
13// Reset only last component which is different
14// FIXME: this is dirty hack, would be good if such functionality will be part
15// of petsc
16
17static PetscErrorCode ierr;
18
19#undef __FUNCT__
20#define __FUNCT__ "only_last_pc_reset_mg"
21PetscErrorCode only_last_pc_reset_mg(PC pc) {
22 PC_MG *mg = (PC_MG *)pc->data;
23 PC_MG_Levels **mglevels = mg->levels;
24
25 PetscInt i, n;
26
27 PetscFunctionBegin;
28 if (mglevels) {
29 n = mglevels[0]->levels;
30
31 for (i = n - 2; i < n - 1; i++) {
32 if (i < 0)
33 continue;
34 CHKERR MatDestroy(&mglevels[i + 1]->restrct);
35 CHKERR MatDestroy(&mglevels[i + 1]->interpolate);
36 CHKERR VecDestroy(&mglevels[i + 1]->rscale);
37 mglevels[i + 1]->restrct = PETSC_NULL;
38 mglevels[i + 1]->interpolate = PETSC_NULL;
39 mglevels[i + 1]->rscale = PETSC_NULL;
40 }
41
42 for (i = n - 2; i < n - 1; i++) {
43 if (i < 0)
44 continue;
45 CHKERR VecDestroy(&mglevels[i + 1]->r);
46 CHKERR VecDestroy(&mglevels[i]->b);
47 CHKERR VecDestroy(&mglevels[i]->x);
48 }
49
50 for (i = 0; i < n; i++) {
51 CHKERR MatDestroy(&mglevels[i]->A);
52 }
53
54 for (i = 0; i < n; i++) {
55 if (i <= 0)
56 continue;
57 // CHKERR MatDestroy(&mglevels[i]->A);CHKERRQ(ierr);
58 if (mglevels[i]->smoothd != mglevels[i]->smoothu) {
59 CHKERR KSPReset(mglevels[i]->smoothd);
60 }
61 CHKERR KSPReset(mglevels[i]->smoothu);
62 }
63 }
64 PetscFunctionReturn(0);
65}
66
67// Create all mg levels or reset last and add one.
68// FIXME: this is dirty hack, would be good if such functionality will be part
69// of petsc
70
71#undef __FUNCT__
72#define __FUNCT__ "pc_mg_set_last_level"
73PetscErrorCode pc_mg_set_last_level(PC pc, PetscInt levels, MPI_Comm *comms) {
74
75 PC_MG *mg = (PC_MG *)pc->data;
76 MPI_Comm comm;
77 PC_MG_Levels **mglevels = mg->levels;
78 PCMGType mgtype = mg->am;
79 PetscInt mgctype = (PetscInt)PC_MG_CYCLE_V;
80 PetscInt i;
81 PetscMPIInt size;
82 const char *prefix;
83 PC ipc;
84 PetscInt n;
85
86 PetscFunctionBegin;
87 PetscValidHeaderSpecific(pc, PC_CLASSID, 1);
88 PetscValidLogicalCollectiveInt(pc, levels, 2);
89 CHKERR PetscObjectGetComm((PetscObject)pc, &comm);
90 vector<bool> set_ksp(levels, false);
91 if (mglevels) {
92 mgctype = mglevels[0]->cycles;
93 /* free up the last level */
95 n = mglevels[0]->levels;
96 if (levels > 2) {
97 set_ksp[0] = true;
98 for (i = 1; i < n - 1; i++) {
99 set_ksp[i] = true;
100 }
101 }
102 for (i = 0; i < n; i++) {
103 if (set_ksp[i])
104 continue;
105 if (mglevels[i]->smoothd != mglevels[i]->smoothu) {
106 CHKERR KSPDestroy(&mglevels[i]->smoothd);
107 }
108 CHKERR KSPDestroy(&mglevels[i]->smoothu);
109 CHKERR PetscFree(mglevels[i]);
110 }
111 }
112
113 mg->nlevels = levels;
114
115 CHKERR PetscMalloc1(levels, &mglevels);
116 CHKERR PetscLogObjectMemory((PetscObject)pc, levels * (sizeof(PC_MG *)));
117
118 CHKERR PCGetOptionsPrefix(pc, &prefix);
119
120 mg->stageApply = 0;
121 for (i = 0; i < levels; i++) {
122 CHKERR PetscNewLog(pc, &mglevels[i]);
123
124 mglevels[i]->levels = levels;
125
126 if (set_ksp[i]) {
127 // Copy from old one to new one
128 if (i != mg->levels[i]->level) {
129 SETERRQ(comm, MOFEM_DATA_INCONSISTENCY, "data inconsistency");
130 }
131 mglevels[i]->level = mg->levels[i]->level;
132 mglevels[i]->cycles = mg->levels[i]->cycles;
133 mglevels[i]->eventsmoothsetup = mg->levels[i]->eventsmoothsetup;
134 mglevels[i]->eventsmoothsolve = mg->levels[i]->eventsmoothsolve;
135 mglevels[i]->eventresidual = mg->levels[i]->eventresidual;
136 mglevels[i]->eventinterprestrict = mg->levels[i]->eventinterprestrict;
137 mglevels[i]->restrct = mg->levels[i]->restrct;
138 mglevels[i]->interpolate = mg->levels[i]->interpolate;
139 mglevels[i]->rscale = mg->levels[i]->rscale;
140 mglevels[i]->r = mg->levels[i]->r;
141 mglevels[i]->b = mg->levels[i]->b;
142 mglevels[i]->x = mg->levels[i]->x;
143 } else {
144 mglevels[i]->level = i;
145 mglevels[i]->cycles = mgctype;
146 mg->default_smoothu = 2;
147 mg->default_smoothd = 2;
148 mglevels[i]->eventsmoothsetup = 0;
149 mglevels[i]->eventsmoothsolve = 0;
150 mglevels[i]->eventresidual = 0;
151 mglevels[i]->eventinterprestrict = 0;
152 }
153
154 if (set_ksp[i]) {
155
156 mglevels[i]->smoothd = mg->levels[i]->smoothd;
157 mglevels[i]->smoothu = mg->levels[i]->smoothu;
158
159 } else {
160
161 // Create all or create last two
162 if (comms)
163 comm = comms[i];
164 CHKERR KSPCreate(comm, &mglevels[i]->smoothd);
165 CHKERR
166 KSPSetErrorIfNotConverged(mglevels[i]->smoothd, pc->erroriffailure);
167 CHKERR KSPSetType(mglevels[i]->smoothd, KSPCHEBYSHEV);
168 CHKERR KSPSetConvergenceTest(mglevels[i]->smoothd, KSPConvergedSkip, NULL,
169 NULL);
170 CHKERR KSPSetNormType(mglevels[i]->smoothd, KSP_NORM_NONE);
171 CHKERR KSPGetPC(mglevels[i]->smoothd, &ipc);
172 CHKERR PCSetType(ipc, PCSOR);
173 CHKERR PetscObjectIncrementTabLevel((PetscObject)mglevels[i]->smoothd,
174 (PetscObject)pc, levels - i);
175 CHKERR
176 KSPSetTolerances(mglevels[i]->smoothd, PETSC_DEFAULT, PETSC_DEFAULT,
177 PETSC_DEFAULT, i ? mg->default_smoothd : 1);
178 CHKERR KSPSetOptionsPrefix(mglevels[i]->smoothd, prefix);
179
180 /* do special stuff for coarse grid */
181 if (!i && levels > 1) {
182 CHKERR KSPAppendOptionsPrefix(mglevels[0]->smoothd, "mg_coarse_");
183
184 /* coarse solve is (redundant) LU by default; set shifttype NONZERO to
185 * avoid annoying zero-pivot in LU preconditioner */
186 CHKERR KSPSetType(mglevels[0]->smoothd, KSPPREONLY);
187 CHKERR KSPGetPC(mglevels[0]->smoothd, &ipc);
188 CHKERR MPI_Comm_size(comm, &size);
189 if (size > 1) {
190 KSP innerksp;
191 PC innerpc;
192 CHKERR PCSetType(ipc, PCREDUNDANT);
193 CHKERR PCRedundantGetKSP(ipc, &innerksp);
194 CHKERR KSPGetPC(innerksp, &innerpc);
195 CHKERR PCFactorSetShiftType(innerpc, MAT_SHIFT_INBLOCKS);
196 } else {
197 CHKERR PCSetType(ipc, PCLU);
198 CHKERR PCFactorSetShiftType(ipc, MAT_SHIFT_INBLOCKS);
199 }
200
201 } else {
202 char tprefix[128];
203 sprintf(tprefix, "mg_levels_%d_", (int)i);
204 CHKERR KSPAppendOptionsPrefix(mglevels[i]->smoothd, tprefix);
205 }
206 CHKERR PetscLogObjectParent((PetscObject)pc,
207 (PetscObject)mglevels[i]->smoothd);
208
209 mglevels[i]->smoothu = mglevels[i]->smoothd;
210 mg->rtol = 0.0;
211 mg->abstol = 0.0;
212 mg->dtol = 0.0;
213 mg->ttol = 0.0;
214 mg->cyclesperpcapply = 1;
215 }
216 }
217 if (mg->levels) {
218 CHKERR PetscFree(mg->levels);
219 }
220 mg->levels = mglevels;
221 CHKERR PCMGSetType(pc, mgtype);
222 PetscFunctionReturn(0);
223}
224
225#undef __FUNCT__
226#define __FUNCT__ "pc_setup_mg"
227PetscErrorCode pc_setup_mg(PC pc) {
228 PC_MG *mg = (PC_MG *)pc->data;
229 PC_MG_Levels **mglevels = mg->levels;
230
231 PetscInt i, n = mglevels[0]->levels;
232 PC cpc;
233 PetscBool preonly, lu, redundant, cholesky, svd,
234 dump = PETSC_FALSE, opsset, use_amat, missinginterpolate = PETSC_FALSE;
235 Vec tvec;
236 DM *dms;
237 PetscViewer viewer = 0;
238
239 PetscFunctionBegin;
240 /* FIX: Move this to PCSetFromOptions_MG? */
241 if (mg->usedmfornumberoflevels) {
242 PetscInt levels;
243 CHKERR DMGetRefineLevel(pc->dm, &levels);
244 levels++;
245 if (levels > n) { /* the problem is now being solved on a finer grid */
246 CHKERR PCMGSetLevels(pc, levels, NULL);
247 n = levels;
248 CHKERR PCSetFromOptions(pc); /* it is bad to call this here, but otherwise
249 will never be called for the new hierarchy */
250 mglevels = mg->levels;
251 }
252 }
253 CHKERR KSPGetPC(mglevels[0]->smoothd, &cpc);
254
255 /* If user did not provide fine grid operators OR operator was not updated
256 * since last global KSPSetOperators() */
257 /* so use those from global PC */
258 /* Is this what we always want? What if user wants to keep old one? */
259 CHKERR KSPGetOperatorsSet(mglevels[n - 1]->smoothd, NULL, &opsset);
260 if (opsset) {
261 Mat mmat;
262 CHKERR KSPGetOperators(mglevels[n - 1]->smoothd, NULL, &mmat);
263 if (mmat == pc->pmat)
264 opsset = PETSC_FALSE;
265 }
266
267 if (!opsset) {
268 CHKERR PCGetUseAmat(pc, &use_amat);
269 if (use_amat) {
270 CHKERR PetscInfo(
271 pc,
272 "Using outer operators to define finest grid operator \n because "
273 "PCMGGetSmoother(pc,nlevels-1,&ksp);KSPSetOperators(ksp,...); was "
274 "not called.\n");
275 CHKERR KSPSetOperators(mglevels[n - 1]->smoothd, pc->mat, pc->pmat);
276 } else {
277 CHKERR PetscInfo(pc, "Using matrix (pmat) operators to define finest "
278 "grid operator \n because "
279 "PCMGGetSmoother(pc,nlevels-1,&ksp);KSPSetOperators("
280 "ksp,...); was not called.\n");
281 CHKERR KSPSetOperators(mglevels[n - 1]->smoothd, pc->pmat, pc->pmat);
282 }
283 }
284
285 for (i = n - 1; i > 0; i--) {
286 if (!(mglevels[i]->interpolate || mglevels[i]->restrct)) {
287 missinginterpolate = PETSC_TRUE;
288 continue;
289 }
290 }
291
292 vector<bool> set_dms_d(n, false);
293 vector<bool> set_dms_u(n, false);
294 for (int i = 0; i < n; i++) {
295 if (mglevels[i]->smoothd->dm)
296 set_dms_d[i] = true;
297 if (mglevels[i]->smoothu->dm)
298 set_dms_u[i] = true;
299 }
300
301 /*
302 Skipping if user has provided all interpolation/restriction needed (since DM
303 might not be able to produce them (when coming from SNES/TS) Skipping for
304 galerkin==2 (externally managed hierarchy such as ML and GAMG). Cleaner logic
305 here would be great. Wrap ML/GAMG as DMs?
306 */
307 if (missinginterpolate && pc->dm && mg->galerkin != 2 && !pc->setupcalled) {
308 /* construct the interpolation from the DMs */
309 Mat p;
310 Vec rscale;
311 CHKERR PetscMalloc1(n, &dms);
312 dms[n - 1] = pc->dm;
313 /* Separately create them so we do not get DMKSP interference between levels
314 */
315 for (i = n - 2; i > -1; i--) {
316 if (!set_dms_d[i]) {
317 CHKERR DMCoarsen(dms[i + 1], MPI_COMM_NULL, &dms[i]);
318 CHKERR KSPSetDM(mglevels[i]->smoothd, dms[i]);
319 if (mg->galerkin) {
320 CHKERR KSPSetDMActive(mglevels[i]->smoothd, PETSC_FALSE);
321 }
322 DMKSP kdm;
323 CHKERR DMGetDMKSPWrite(dms[i], &kdm);
324 /* Ugly hack so that the next KSPSetUp() will use the RHS that we set. A
325 * better fix is to change dmActive to take a bitwise OR of computing
326 * the matrix, RHS, and initial iterate. */
327 kdm->ops->computerhs = NULL;
328 kdm->rhsctx = NULL;
329 } else {
330 dms[i] = mglevels[i]->smoothd->dm;
331 dms[i]->leveldown = dms[i + 1]->leveldown + 1;
332 CHKERR PetscObjectReference((PetscObject)dms[i]);
333 }
334 }
335 for (i = n - 2; i > -1; i--) {
336 if (!mglevels[i + 1]->interpolate) {
337 // cerr << "interpolate " << i+1 << endl;
338 CHKERR DMCreateInterpolation(dms[i], dms[i + 1], &p, &rscale);
339 CHKERR PCMGSetInterpolation(pc, i + 1, p);
340 if (rscale) {
341 CHKERR PCMGSetRScale(pc, i + 1, rscale);
342 }
343 CHKERR VecDestroy(&rscale);
344 CHKERR MatDestroy(&p);
345 } else {
346 // cerr << "not interpolation " << i+1 << endl;
347 }
348 }
349 for (i = n - 2; i > -1; i--) {
350 CHKERR DMDestroy(&dms[i]);
351 }
352 CHKERR PetscFree(dms);
353 }
354
355 if (pc->dm && !pc->setupcalled) {
356 /* finest smoother also gets DM but it is not active, independent of whether
357 * galerkin==2 */
358 CHKERR KSPSetDM(mglevels[n - 1]->smoothd, pc->dm);
359 CHKERR KSPSetDMActive(mglevels[n - 1]->smoothd, PETSC_FALSE);
360 }
361
362 if (mg->galerkin == 1) {
363 MPI_Comm comm;
364 CHKERR PetscObjectGetComm((PetscObject)pc, &comm);
365 SETERRQ(comm, MOFEM_DATA_INCONSISTENCY, "should not happen");
366 } // else if (!mg->galerkin && pc->dm && pc->dm->x) {
367 // /* need to restrict Jacobian location to coarser meshes for evaluation */
368 // for (i=n-2; i>-1; i--) {
369 // Mat R;
370 // Vec rscale;
371 // if (!mglevels[i]->smoothd->dm->x) {
372 // Vec *vecs;
373 // CHKERR
374 // KSPCreateVecs(mglevels[i]->smoothd,1,&vecs,0,NULL);CHKERRQ(ierr);
375 //
376 // mglevels[i]->smoothd->dm->x = vecs[0];
377 //
378 // CHKERR PetscFree(vecs);CHKERRQ(ierr);
379 // }
380 // CHKERR PCMGGetRestriction(pc,i+1,&R);CHKERRQ(ierr);
381 // CHKERR PCMGGetRScale(pc,i+1,&rscale);CHKERRQ(ierr);
382 // CHKERR
383 // MatRestrict(R,mglevels[i+1]->smoothd->dm->x,mglevels[i]->smoothd->dm->x);CHKERRQ(ierr);
384 // CHKERR
385 // VecPointwiseMult(mglevels[i]->smoothd->dm->x,mglevels[i]->smoothd->dm->x,rscale);CHKERRQ(ierr);
386 // }
387 // }
388 if (!mg->galerkin && pc->dm) {
389 for (i = n - 2; i >= 0; i--) {
390 if (!mglevels[i + 1]->restrct) {
391 // cerr << n << " " << i << endl;
392 DM dmfine, dmcoarse;
393 Mat Restrict, Inject;
394 Vec rscale;
395 CHKERR KSPGetDM(mglevels[i + 1]->smoothd, &dmfine);
396 CHKERR KSPGetDM(mglevels[i]->smoothd, &dmcoarse);
397 CHKERR PCMGGetRestriction(pc, i + 1, &Restrict);
398 CHKERR PCMGGetRScale(pc, i + 1, &rscale);
399 Inject = NULL; /* Callback should create it if it needs Injection */
400 CHKERR DMRestrict(dmfine, Restrict, rscale, Inject, dmcoarse);
401 }
402 }
403 }
404
405 if (!pc->setupcalled) {
406 for (i = 0; i < n; i++) {
407 if (set_dms_d[i])
408 continue;
409 CHKERR KSPSetFromOptions(mglevels[i]->smoothd);
410 }
411 for (i = 1; i < n; i++) {
412 if (set_dms_u[i])
413 continue;
414 if (mglevels[i]->smoothu &&
415 (mglevels[i]->smoothu != mglevels[i]->smoothd)) {
416 CHKERR KSPSetFromOptions(mglevels[i]->smoothu);
417 }
418 }
419 for (i = 1; i < n; i++) {
420 CHKERR PCMGGetInterpolation(pc, i, &mglevels[i]->interpolate);
421 CHKERR PCMGGetRestriction(pc, i, &mglevels[i]->restrct);
422 }
423 for (i = 0; i < n - 1; i++) {
424 if (!mglevels[i]->b) {
425 Vec *vec;
426 CHKERR KSPCreateVecs(mglevels[i]->smoothd, 1, &vec, 0, NULL);
427 CHKERR PCMGSetRhs(pc, i, *vec);
428 CHKERR VecDestroy(vec);
429 CHKERR PetscFree(vec);
430 }
431 if (!mglevels[i]->r && i) {
432 CHKERR VecDuplicate(mglevels[i]->b, &tvec);
433 CHKERR PCMGSetR(pc, i, tvec);
434 CHKERR VecDestroy(&tvec);
435 }
436 if (!mglevels[i]->x) {
437 CHKERR VecDuplicate(mglevels[i]->b, &tvec);
438 CHKERR PCMGSetX(pc, i, tvec);
439 CHKERR VecDestroy(&tvec);
440 }
441 }
442 if (n != 1 && !mglevels[n - 1]->r) {
443 /* PCMGSetR() on the finest level if user did not supply it */
444 Vec *vec;
445 CHKERR KSPCreateVecs(mglevels[n - 1]->smoothd, 1, &vec, 0, NULL);
446 CHKERR PCMGSetR(pc, n - 1, *vec);
447 CHKERR VecDestroy(vec);
448 CHKERR PetscFree(vec);
449 }
450 }
451
452 if (pc->dm) {
453 /* need to tell all the coarser levels to rebuild the matrix using the DM
454 * for that level */
455 for (i = 0; i < n - 1; i++) {
456 if (set_dms_d[i])
457 continue;
458 if (mglevels[i]->smoothd->setupstage != KSP_SETUP_NEW)
459 mglevels[i]->smoothd->setupstage = KSP_SETUP_NEWMATRIX;
460 }
461 }
462
463 for (i = 1; i < n; i++) {
464 if (mglevels[i]->smoothu == mglevels[i]->smoothd || mg->am == PC_MG_FULL ||
465 mg->am == PC_MG_KASKADE || mg->cyclesperpcapply > 1) {
466 /* if doing only down then initial guess is zero */
467 CHKERR KSPSetInitialGuessNonzero(mglevels[i]->smoothd, PETSC_TRUE);
468 }
469 if (mglevels[i]->eventsmoothsetup) {
470 CHKERR PetscLogEventBegin(mglevels[i]->eventsmoothsetup, 0, 0, 0, 0);
471 }
472 CHKERR KSPSetUp(mglevels[i]->smoothd);
473 if (mglevels[i]->eventsmoothsetup) {
474 CHKERR PetscLogEventEnd(mglevels[i]->eventsmoothsetup, 0, 0, 0, 0);
475 }
476 if (!mglevels[i]->residual) {
477 Mat mat;
478 CHKERR KSPGetOperators(mglevels[i]->smoothd, NULL, &mat);
479 CHKERR PCMGSetResidual(pc, i, PCMGResidualDefault, mat);
480 }
481 }
482 for (i = 1; i < n; i++) {
483 if (mglevels[i]->smoothu && mglevels[i]->smoothu != mglevels[i]->smoothd) {
484 Mat downmat, downpmat;
485
486 /* check if operators have been set for up, if not use down operators to
487 * set them */
488 if (!set_dms_u[i]) {
489 CHKERR KSPGetOperatorsSet(mglevels[i]->smoothu, &opsset, NULL);
490 if (!opsset) {
491 CHKERR KSPGetOperators(mglevels[i]->smoothd, &downmat, &downpmat);
492 CHKERR KSPSetOperators(mglevels[i]->smoothu, downmat, downpmat);
493 }
494 CHKERR KSPSetInitialGuessNonzero(mglevels[i]->smoothu, PETSC_TRUE);
495 if (mglevels[i]->eventsmoothsetup) {
496 CHKERR PetscLogEventBegin(mglevels[i]->eventsmoothsetup, 0, 0, 0, 0);
497 }
498 CHKERR KSPSetUp(mglevels[i]->smoothu);
499 if (mglevels[i]->eventsmoothsetup) {
500 CHKERR PetscLogEventEnd(mglevels[i]->eventsmoothsetup, 0, 0, 0, 0);
501 }
502 }
503 }
504 }
505
506 /*
507 If coarse solver is not direct method then DO NOT USE preonly
508 */
509 CHKERR PetscObjectTypeCompare((PetscObject)mglevels[0]->smoothd, KSPPREONLY,
510 &preonly);
511 if (preonly) {
512 CHKERR PetscObjectTypeCompare((PetscObject)cpc, PCLU, &lu);
513 CHKERR PetscObjectTypeCompare((PetscObject)cpc, PCREDUNDANT, &redundant);
514 CHKERR PetscObjectTypeCompare((PetscObject)cpc, PCCHOLESKY, &cholesky);
515 CHKERR PetscObjectTypeCompare((PetscObject)cpc, PCSVD, &svd);
516 if (!lu && !redundant && !cholesky && !svd) {
517 CHKERR KSPSetType(mglevels[0]->smoothd, KSPGMRES);
518 }
519 }
520
521 if (mglevels[0]->eventsmoothsetup) {
522 CHKERR PetscLogEventBegin(mglevels[0]->eventsmoothsetup, 0, 0, 0, 0);
523 }
524 CHKERR KSPSetUp(mglevels[0]->smoothd);
525 if (mglevels[0]->eventsmoothsetup) {
526 CHKERR PetscLogEventEnd(mglevels[0]->eventsmoothsetup, 0, 0, 0, 0);
527 }
528
529 /*
530 Dump the interpolation/restriction matrices plus the
531 Jacobian/stiffness on each level. This allows MATLAB users to
532 easily check if the Galerkin condition A_c = R A_f R^T is satisfied.
533
534 Only support one or the other at the same time.
535 */
536 // #if defined(PETSC_USE_SOCKET_VIEWER)
537 // CHKERR
538 // PetscOptionsGetBool(((PetscObject)pc)->prefix,"-pc_mg_dump_matlab",&dump,NULL);CHKERRQ(ierr);
539 // if (dump) viewer =
540 // PETSC_VIEWER_SOCKET_(PetscObjectComm((PetscObject)pc)); dump =
541 // PETSC_FALSE;
542 // #endif
543 // CHKERR
544 // PetscOptionsGetBool(((PetscObject)pc)->prefix,"-pc_mg_dump_binary",&dump,NULL);CHKERRQ(ierr);
545 // if (dump) viewer =
546 // PETSC_VIEWER_BINARY_(PetscObjectComm((PetscObject)pc));
547
548 if (viewer) {
549 for (i = 1; i < n; i++) {
550 CHKERR MatView(mglevels[i]->restrct, viewer);
551 }
552 for (i = 0; i < n; i++) {
553 CHKERR KSPGetPC(mglevels[i]->smoothd, &pc);
554 CHKERR MatView(pc->mat, viewer);
555 }
556 }
557 PetscFunctionReturn(0);
558}
559
560} // namespace SolidShellModule
static Index< 'p', 3 > p
@ MOFEM_DATA_INCONSISTENCY
Definition definitions.h:31
#define CHKERR
Inline error check.
FTensor::Index< 'n', SPACE_DIM > n
FTensor::Index< 'i', SPACE_DIM > i
PetscErrorCode pc_setup_mg(PC pc)
Definition HackMG.cpp:227
PetscErrorCode pc_mg_set_last_level(PC pc, PetscInt levels, MPI_Comm *comms)
Definition HackMG.cpp:73
PetscErrorCode only_last_pc_reset_mg(PC pc)
Definition HackMG.cpp:21
static PetscErrorCode ierr
Definition HackMG.cpp:17