Actual source code: pcis.c


  2: #include <petsc/private/pcisimpl.h>

  4: static PetscErrorCode PCISSetUseStiffnessScaling_IS(PC pc, PetscBool use)
  5: {
  6:   PC_IS *pcis = (PC_IS *)pc->data;

  8:   pcis->use_stiffness_scaling = use;
  9:   return 0;
 10: }

 12: /*@
 13:    PCISSetUseStiffnessScaling - Tells `PCIS` to construct partition of unity using
 14:                               the local matrices' diagonal entries

 16:    Logically Collective

 18:    Input Parameters:
 19: +  pc - the preconditioning context
 20: -  use - whether or not pcis use matrix diagonal to build partition of unity.

 22:    Level: intermediate

 24:    Developer Note:
 25:    There is no manual page for `PCIS` nor some of its methods

 27: .seealso: `PCIS`, `PCBDDC`
 28: @*/
 29: PetscErrorCode PCISSetUseStiffnessScaling(PC pc, PetscBool use)
 30: {
 33:   PetscTryMethod(pc, "PCISSetUseStiffnessScaling_C", (PC, PetscBool), (pc, use));
 34:   return 0;
 35: }

 37: static PetscErrorCode PCISSetSubdomainDiagonalScaling_IS(PC pc, Vec scaling_factors)
 38: {
 39:   PC_IS *pcis = (PC_IS *)pc->data;

 41:   PetscObjectReference((PetscObject)scaling_factors);
 42:   VecDestroy(&pcis->D);
 43:   pcis->D = scaling_factors;
 44:   if (pc->setupcalled) {
 45:     PetscInt sn;

 47:     VecGetSize(pcis->D, &sn);
 48:     if (sn == pcis->n) {
 49:       VecScatterBegin(pcis->N_to_B, pcis->D, pcis->vec1_B, INSERT_VALUES, SCATTER_FORWARD);
 50:       VecScatterEnd(pcis->N_to_B, pcis->D, pcis->vec1_B, INSERT_VALUES, SCATTER_FORWARD);
 51:       VecDestroy(&pcis->D);
 52:       VecDuplicate(pcis->vec1_B, &pcis->D);
 53:       VecCopy(pcis->vec1_B, pcis->D);
 55:   }
 56:   return 0;
 57: }

 59: /*@
 60:    PCISSetSubdomainDiagonalScaling - Set diagonal scaling for `PCIS`.

 62:    Logically Collective

 64:    Input Parameters:
 65: +  pc - the preconditioning context
 66: -  scaling_factors - scaling factors for the subdomain

 68:    Level: intermediate

 70:    Note:
 71:    Intended for use with jumping coefficients cases.

 73:    Developer Note:
 74:    There is no manual page for `PCIS` nor some of its methods

 76: .seealso: `PCIS`, `PCBDDC`
 77: @*/
 78: PetscErrorCode PCISSetSubdomainDiagonalScaling(PC pc, Vec scaling_factors)
 79: {
 82:   PetscTryMethod(pc, "PCISSetSubdomainDiagonalScaling_C", (PC, Vec), (pc, scaling_factors));
 83:   return 0;
 84: }

 86: static PetscErrorCode PCISSetSubdomainScalingFactor_IS(PC pc, PetscScalar scal)
 87: {
 88:   PC_IS *pcis = (PC_IS *)pc->data;

 90:   pcis->scaling_factor = scal;
 91:   if (pcis->D) VecSet(pcis->D, pcis->scaling_factor);
 92:   return 0;
 93: }

 95: /*@
 96:  PCISSetSubdomainScalingFactor - Set scaling factor for `PCIS`.

 98:    Not collective

100:    Input Parameters:
101: +  pc - the preconditioning context
102: -  scal - scaling factor for the subdomain

104:    Level: intermediate

106:    Note:
107:    Intended for use with the jumping coefficients cases.

109:    Developer Note:
110:    There is no manual page for `PCIS` nor some of its methods

112: .seealso: `PCIS`, `PCBDDC`
113: @*/
114: PetscErrorCode PCISSetSubdomainScalingFactor(PC pc, PetscScalar scal)
115: {
117:   PetscTryMethod(pc, "PCISSetSubdomainScalingFactor_C", (PC, PetscScalar), (pc, scal));
118:   return 0;
119: }

121: PetscErrorCode PCISSetUp(PC pc, PetscBool computematrices, PetscBool computesolvers)
122: {
123:   PC_IS    *pcis = (PC_IS *)(pc->data);
124:   Mat_IS   *matis;
125:   MatReuse  reuse;
126:   PetscBool flg, issbaij;

128:   PetscObjectTypeCompare((PetscObject)pc->pmat, MATIS, &flg);
130:   matis = (Mat_IS *)pc->pmat->data;
131:   if (pc->useAmat) {
132:     PetscObjectTypeCompare((PetscObject)pc->mat, MATIS, &flg);
134:   }

136:   /* first time creation, get info on substructuring */
137:   if (!pc->setupcalled) {
138:     PetscInt  n_I;
139:     PetscInt *idx_I_local, *idx_B_local, *idx_I_global, *idx_B_global;
140:     PetscBT   bt;
141:     PetscInt  i, j;

143:     /* get info on mapping */
144:     PetscObjectReference((PetscObject)matis->rmapping);
145:     ISLocalToGlobalMappingDestroy(&pcis->mapping);
146:     pcis->mapping = matis->rmapping;
147:     ISLocalToGlobalMappingGetSize(pcis->mapping, &pcis->n);
148:     ISLocalToGlobalMappingGetInfo(pcis->mapping, &(pcis->n_neigh), &(pcis->neigh), &(pcis->n_shared), &(pcis->shared));

150:     /* Identifying interior and interface nodes, in local numbering */
151:     PetscBTCreate(pcis->n, &bt);
152:     for (i = 0; i < pcis->n_neigh; i++)
153:       for (j = 0; j < pcis->n_shared[i]; j++) PetscBTSet(bt, pcis->shared[i][j]);

155:     /* Creating local and global index sets for interior and interface nodes. */
156:     PetscMalloc1(pcis->n, &idx_I_local);
157:     PetscMalloc1(pcis->n, &idx_B_local);
158:     for (i = 0, pcis->n_B = 0, n_I = 0; i < pcis->n; i++) {
159:       if (!PetscBTLookup(bt, i)) {
160:         idx_I_local[n_I] = i;
161:         n_I++;
162:       } else {
163:         idx_B_local[pcis->n_B] = i;
164:         pcis->n_B++;
165:       }
166:     }

168:     /* Getting the global numbering */
169:     idx_B_global = idx_I_local + n_I; /* Just avoiding allocating extra memory, since we have vacant space */
170:     idx_I_global = idx_B_local + pcis->n_B;
171:     ISLocalToGlobalMappingApply(pcis->mapping, pcis->n_B, idx_B_local, idx_B_global);
172:     ISLocalToGlobalMappingApply(pcis->mapping, n_I, idx_I_local, idx_I_global);

174:     /* Creating the index sets */
175:     ISCreateGeneral(PETSC_COMM_SELF, pcis->n_B, idx_B_local, PETSC_COPY_VALUES, &pcis->is_B_local);
176:     ISCreateGeneral(PetscObjectComm((PetscObject)pc), pcis->n_B, idx_B_global, PETSC_COPY_VALUES, &pcis->is_B_global);
177:     ISCreateGeneral(PETSC_COMM_SELF, n_I, idx_I_local, PETSC_COPY_VALUES, &pcis->is_I_local);
178:     ISCreateGeneral(PetscObjectComm((PetscObject)pc), n_I, idx_I_global, PETSC_COPY_VALUES, &pcis->is_I_global);

180:     /* Freeing memory */
181:     PetscFree(idx_B_local);
182:     PetscFree(idx_I_local);
183:     PetscBTDestroy(&bt);

185:     /* Creating work vectors and arrays */
186:     VecDuplicate(matis->x, &pcis->vec1_N);
187:     VecDuplicate(pcis->vec1_N, &pcis->vec2_N);
188:     VecCreate(PETSC_COMM_SELF, &pcis->vec1_D);
189:     VecSetSizes(pcis->vec1_D, pcis->n - pcis->n_B, PETSC_DECIDE);
190:     VecSetType(pcis->vec1_D, ((PetscObject)pcis->vec1_N)->type_name);
191:     VecDuplicate(pcis->vec1_D, &pcis->vec2_D);
192:     VecDuplicate(pcis->vec1_D, &pcis->vec3_D);
193:     VecDuplicate(pcis->vec1_D, &pcis->vec4_D);
194:     VecCreate(PETSC_COMM_SELF, &pcis->vec1_B);
195:     VecSetSizes(pcis->vec1_B, pcis->n_B, PETSC_DECIDE);
196:     VecSetType(pcis->vec1_B, ((PetscObject)pcis->vec1_N)->type_name);
197:     VecDuplicate(pcis->vec1_B, &pcis->vec2_B);
198:     VecDuplicate(pcis->vec1_B, &pcis->vec3_B);
199:     MatCreateVecs(pc->pmat, &pcis->vec1_global, NULL);
200:     PetscMalloc1(pcis->n, &pcis->work_N);
201:     /* scaling vector */
202:     if (!pcis->D) { /* it can happen that the user passed in a scaling vector via PCISSetSubdomainDiagonalScaling */
203:       VecDuplicate(pcis->vec1_B, &pcis->D);
204:       VecSet(pcis->D, pcis->scaling_factor);
205:     }

207:     /* Creating the scatter contexts */
208:     VecScatterCreate(pcis->vec1_N, pcis->is_I_local, pcis->vec1_D, (IS)0, &pcis->N_to_D);
209:     VecScatterCreate(pcis->vec1_global, pcis->is_I_global, pcis->vec1_D, (IS)0, &pcis->global_to_D);
210:     VecScatterCreate(pcis->vec1_N, pcis->is_B_local, pcis->vec1_B, (IS)0, &pcis->N_to_B);
211:     VecScatterCreate(pcis->vec1_global, pcis->is_B_global, pcis->vec1_B, (IS)0, &pcis->global_to_B);

213:     /* map from boundary to local */
214:     ISLocalToGlobalMappingCreateIS(pcis->is_B_local, &pcis->BtoNmap);
215:   }

217:   {
218:     PetscInt sn;

220:     VecGetSize(pcis->D, &sn);
221:     if (sn == pcis->n) {
222:       VecScatterBegin(pcis->N_to_B, pcis->D, pcis->vec1_B, INSERT_VALUES, SCATTER_FORWARD);
223:       VecScatterEnd(pcis->N_to_B, pcis->D, pcis->vec1_B, INSERT_VALUES, SCATTER_FORWARD);
224:       VecDestroy(&pcis->D);
225:       VecDuplicate(pcis->vec1_B, &pcis->D);
226:       VecCopy(pcis->vec1_B, pcis->D);
228:   }

230:   /*
231:     Extracting the blocks A_II, A_BI, A_IB and A_BB from A. If the numbering
232:     is such that interior nodes come first than the interface ones, we have

234:         [ A_II | A_IB ]
235:     A = [------+------]
236:         [ A_BI | A_BB ]
237:   */
238:   if (computematrices) {
239:     PetscBool amat = (PetscBool)(pc->mat != pc->pmat && pc->useAmat);
240:     PetscInt  bs, ibs;

242:     reuse = MAT_INITIAL_MATRIX;
243:     if (pcis->reusesubmatrices && pc->setupcalled) {
244:       if (pc->flag == SAME_NONZERO_PATTERN) {
245:         reuse = MAT_REUSE_MATRIX;
246:       } else {
247:         reuse = MAT_INITIAL_MATRIX;
248:       }
249:     }
250:     if (reuse == MAT_INITIAL_MATRIX) {
251:       MatDestroy(&pcis->A_II);
252:       MatDestroy(&pcis->pA_II);
253:       MatDestroy(&pcis->A_IB);
254:       MatDestroy(&pcis->A_BI);
255:       MatDestroy(&pcis->A_BB);
256:     }

258:     ISLocalToGlobalMappingGetBlockSize(pcis->mapping, &ibs);
259:     MatGetBlockSize(matis->A, &bs);
260:     MatCreateSubMatrix(matis->A, pcis->is_I_local, pcis->is_I_local, reuse, &pcis->pA_II);
261:     if (amat) {
262:       Mat_IS *amatis = (Mat_IS *)pc->mat->data;
263:       MatCreateSubMatrix(amatis->A, pcis->is_I_local, pcis->is_I_local, reuse, &pcis->A_II);
264:     } else {
265:       PetscObjectReference((PetscObject)pcis->pA_II);
266:       MatDestroy(&pcis->A_II);
267:       pcis->A_II = pcis->pA_II;
268:     }
269:     MatSetBlockSize(pcis->A_II, bs == ibs ? bs : 1);
270:     MatSetBlockSize(pcis->pA_II, bs == ibs ? bs : 1);
271:     MatCreateSubMatrix(matis->A, pcis->is_B_local, pcis->is_B_local, reuse, &pcis->A_BB);
272:     PetscObjectTypeCompare((PetscObject)matis->A, MATSEQSBAIJ, &issbaij);
273:     if (!issbaij) {
274:       MatCreateSubMatrix(matis->A, pcis->is_I_local, pcis->is_B_local, reuse, &pcis->A_IB);
275:       MatCreateSubMatrix(matis->A, pcis->is_B_local, pcis->is_I_local, reuse, &pcis->A_BI);
276:     } else {
277:       Mat newmat;

279:       MatConvert(matis->A, MATSEQBAIJ, MAT_INITIAL_MATRIX, &newmat);
280:       MatCreateSubMatrix(newmat, pcis->is_I_local, pcis->is_B_local, reuse, &pcis->A_IB);
281:       MatCreateSubMatrix(newmat, pcis->is_B_local, pcis->is_I_local, reuse, &pcis->A_BI);
282:       MatDestroy(&newmat);
283:     }
284:     MatSetBlockSize(pcis->A_BB, bs == ibs ? bs : 1);
285:   }

287:   /* Creating scaling vector D */
288:   PetscOptionsGetBool(((PetscObject)pc)->options, ((PetscObject)pc)->prefix, "-pc_is_use_stiffness_scaling", &pcis->use_stiffness_scaling, NULL);
289:   if (pcis->use_stiffness_scaling) {
290:     PetscScalar *a;
291:     PetscInt     i, n;

293:     if (pcis->A_BB) {
294:       MatGetDiagonal(pcis->A_BB, pcis->D);
295:     } else {
296:       MatGetDiagonal(matis->A, pcis->vec1_N);
297:       VecScatterBegin(pcis->N_to_B, pcis->vec1_N, pcis->D, INSERT_VALUES, SCATTER_FORWARD);
298:       VecScatterEnd(pcis->N_to_B, pcis->vec1_N, pcis->D, INSERT_VALUES, SCATTER_FORWARD);
299:     }
300:     VecAbs(pcis->D);
301:     VecGetLocalSize(pcis->D, &n);
302:     VecGetArray(pcis->D, &a);
303:     for (i = 0; i < n; i++)
304:       if (PetscAbsScalar(a[i]) < PETSC_SMALL) a[i] = 1.0;
305:     VecRestoreArray(pcis->D, &a);
306:   }
307:   VecSet(pcis->vec1_global, 0.0);
308:   VecScatterBegin(pcis->global_to_B, pcis->D, pcis->vec1_global, ADD_VALUES, SCATTER_REVERSE);
309:   VecScatterEnd(pcis->global_to_B, pcis->D, pcis->vec1_global, ADD_VALUES, SCATTER_REVERSE);
310:   VecScatterBegin(pcis->global_to_B, pcis->vec1_global, pcis->vec1_B, INSERT_VALUES, SCATTER_FORWARD);
311:   VecScatterEnd(pcis->global_to_B, pcis->vec1_global, pcis->vec1_B, INSERT_VALUES, SCATTER_FORWARD);
312:   VecPointwiseDivide(pcis->D, pcis->D, pcis->vec1_B);
313:   /* See historical note 01, at the bottom of this file. */

315:   /* Creating the KSP contexts for the local Dirichlet and Neumann problems */
316:   if (computesolvers) {
317:     PC pc_ctx;

319:     pcis->pure_neumann = matis->pure_neumann;
320:     /* Dirichlet */
321:     KSPCreate(PETSC_COMM_SELF, &pcis->ksp_D);
322:     KSPSetErrorIfNotConverged(pcis->ksp_D, pc->erroriffailure);
323:     PetscObjectIncrementTabLevel((PetscObject)pcis->ksp_D, (PetscObject)pc, 1);
324:     KSPSetOperators(pcis->ksp_D, pcis->A_II, pcis->A_II);
325:     KSPSetOptionsPrefix(pcis->ksp_D, "is_localD_");
326:     KSPGetPC(pcis->ksp_D, &pc_ctx);
327:     PCSetType(pc_ctx, PCLU);
328:     KSPSetType(pcis->ksp_D, KSPPREONLY);
329:     KSPSetFromOptions(pcis->ksp_D);
330:     /* the vectors in the following line are dummy arguments, just telling the KSP the vector size. Values are not used */
331:     KSPSetUp(pcis->ksp_D);
332:     /* Neumann */
333:     KSPCreate(PETSC_COMM_SELF, &pcis->ksp_N);
334:     KSPSetErrorIfNotConverged(pcis->ksp_N, pc->erroriffailure);
335:     PetscObjectIncrementTabLevel((PetscObject)pcis->ksp_N, (PetscObject)pc, 1);
336:     KSPSetOperators(pcis->ksp_N, matis->A, matis->A);
337:     KSPSetOptionsPrefix(pcis->ksp_N, "is_localN_");
338:     KSPGetPC(pcis->ksp_N, &pc_ctx);
339:     PCSetType(pc_ctx, PCLU);
340:     KSPSetType(pcis->ksp_N, KSPPREONLY);
341:     KSPSetFromOptions(pcis->ksp_N);
342:     {
343:       PetscBool damp_fixed = PETSC_FALSE, remove_nullspace_fixed = PETSC_FALSE, set_damping_factor_floating = PETSC_FALSE, not_damp_floating = PETSC_FALSE, not_remove_nullspace_floating = PETSC_FALSE;
344:       PetscReal fixed_factor, floating_factor;

346:       PetscOptionsGetReal(((PetscObject)pc_ctx)->options, ((PetscObject)pc_ctx)->prefix, "-pc_is_damp_fixed", &fixed_factor, &damp_fixed);
347:       if (!damp_fixed) fixed_factor = 0.0;
348:       PetscOptionsGetBool(((PetscObject)pc_ctx)->options, ((PetscObject)pc_ctx)->prefix, "-pc_is_damp_fixed", &damp_fixed, NULL);

350:       PetscOptionsGetBool(((PetscObject)pc_ctx)->options, ((PetscObject)pc_ctx)->prefix, "-pc_is_remove_nullspace_fixed", &remove_nullspace_fixed, NULL);

352:       PetscOptionsGetReal(((PetscObject)pc_ctx)->options, ((PetscObject)pc_ctx)->prefix, "-pc_is_set_damping_factor_floating", &floating_factor, &set_damping_factor_floating);
353:       if (!set_damping_factor_floating) floating_factor = 0.0;
354:       PetscOptionsGetBool(((PetscObject)pc_ctx)->options, ((PetscObject)pc_ctx)->prefix, "-pc_is_set_damping_factor_floating", &set_damping_factor_floating, NULL);
355:       if (!set_damping_factor_floating) floating_factor = 1.e-12;

357:       PetscOptionsGetBool(((PetscObject)pc_ctx)->options, ((PetscObject)pc_ctx)->prefix, "-pc_is_not_damp_floating", &not_damp_floating, NULL);

359:       PetscOptionsGetBool(((PetscObject)pc_ctx)->options, ((PetscObject)pc_ctx)->prefix, "-pc_is_not_remove_nullspace_floating", &not_remove_nullspace_floating, NULL);

361:       if (pcis->pure_neumann) { /* floating subdomain */
362:         if (!(not_damp_floating)) {
363:           PCFactorSetShiftType(pc_ctx, MAT_SHIFT_NONZERO);
364:           PCFactorSetShiftAmount(pc_ctx, floating_factor);
365:         }
366:         if (!(not_remove_nullspace_floating)) {
367:           MatNullSpace nullsp;
368:           MatNullSpaceCreate(PETSC_COMM_SELF, PETSC_TRUE, 0, NULL, &nullsp);
369:           MatSetNullSpace(matis->A, nullsp);
370:           MatNullSpaceDestroy(&nullsp);
371:         }
372:       } else { /* fixed subdomain */
373:         if (damp_fixed) {
374:           PCFactorSetShiftType(pc_ctx, MAT_SHIFT_NONZERO);
375:           PCFactorSetShiftAmount(pc_ctx, floating_factor);
376:         }
377:         if (remove_nullspace_fixed) {
378:           MatNullSpace nullsp;
379:           MatNullSpaceCreate(PETSC_COMM_SELF, PETSC_TRUE, 0, NULL, &nullsp);
380:           MatSetNullSpace(matis->A, nullsp);
381:           MatNullSpaceDestroy(&nullsp);
382:         }
383:       }
384:     }
385:     /* the vectors in the following line are dummy arguments, just telling the KSP the vector size. Values are not used */
386:     KSPSetUp(pcis->ksp_N);
387:   }
388:   return 0;
389: }

391: /*
392:    PCISDestroy -
393: */
394: PetscErrorCode PCISDestroy(PC pc)
395: {
396:   PC_IS *pcis;

398:   if (!pc) return 0;
399:   pcis = (PC_IS *)(pc->data);
400:   ISDestroy(&pcis->is_B_local);
401:   ISDestroy(&pcis->is_I_local);
402:   ISDestroy(&pcis->is_B_global);
403:   ISDestroy(&pcis->is_I_global);
404:   MatDestroy(&pcis->A_II);
405:   MatDestroy(&pcis->pA_II);
406:   MatDestroy(&pcis->A_IB);
407:   MatDestroy(&pcis->A_BI);
408:   MatDestroy(&pcis->A_BB);
409:   VecDestroy(&pcis->D);
410:   KSPDestroy(&pcis->ksp_N);
411:   KSPDestroy(&pcis->ksp_D);
412:   VecDestroy(&pcis->vec1_N);
413:   VecDestroy(&pcis->vec2_N);
414:   VecDestroy(&pcis->vec1_D);
415:   VecDestroy(&pcis->vec2_D);
416:   VecDestroy(&pcis->vec3_D);
417:   VecDestroy(&pcis->vec4_D);
418:   VecDestroy(&pcis->vec1_B);
419:   VecDestroy(&pcis->vec2_B);
420:   VecDestroy(&pcis->vec3_B);
421:   VecDestroy(&pcis->vec1_global);
422:   VecScatterDestroy(&pcis->global_to_D);
423:   VecScatterDestroy(&pcis->N_to_B);
424:   VecScatterDestroy(&pcis->N_to_D);
425:   VecScatterDestroy(&pcis->global_to_B);
426:   PetscFree(pcis->work_N);
427:   if (pcis->n_neigh > -1) ISLocalToGlobalMappingRestoreInfo(pcis->mapping, &(pcis->n_neigh), &(pcis->neigh), &(pcis->n_shared), &(pcis->shared));
428:   ISLocalToGlobalMappingDestroy(&pcis->mapping);
429:   ISLocalToGlobalMappingDestroy(&pcis->BtoNmap);
430:   PetscObjectComposeFunction((PetscObject)pc, "PCISSetUseStiffnessScaling_C", NULL);
431:   PetscObjectComposeFunction((PetscObject)pc, "PCISSetSubdomainScalingFactor_C", NULL);
432:   PetscObjectComposeFunction((PetscObject)pc, "PCISSetSubdomainDiagonalScaling_C", NULL);
433:   return 0;
434: }

436: /*
437:    PCISCreate -
438: */
439: PetscErrorCode PCISCreate(PC pc)
440: {
441:   PC_IS *pcis = (PC_IS *)(pc->data);

443:   if (!pcis) {
444:     PetscNew(&pcis);
445:     pc->data = pcis;
446:   }
447:   pcis->n_neigh          = -1;
448:   pcis->scaling_factor   = 1.0;
449:   pcis->reusesubmatrices = PETSC_TRUE;
450:   /* composing functions */
451:   PetscObjectComposeFunction((PetscObject)pc, "PCISSetUseStiffnessScaling_C", PCISSetUseStiffnessScaling_IS);
452:   PetscObjectComposeFunction((PetscObject)pc, "PCISSetSubdomainScalingFactor_C", PCISSetSubdomainScalingFactor_IS);
453:   PetscObjectComposeFunction((PetscObject)pc, "PCISSetSubdomainDiagonalScaling_C", PCISSetSubdomainDiagonalScaling_IS);
454:   return 0;
455: }

457: /*
458:    PCISApplySchur -

460:    Input parameters:
461: .  pc - preconditioner context
462: .  v - vector to which the Schur complement is to be applied (it is NOT modified inside this function, UNLESS vec2_B is null)

464:    Output parameters:
465: .  vec1_B - result of Schur complement applied to chunk
466: .  vec2_B - garbage (used as work space), or null (and v is used as workspace)
467: .  vec1_D - garbage (used as work space)
468: .  vec2_D - garbage (used as work space)

470: */
471: PetscErrorCode PCISApplySchur(PC pc, Vec v, Vec vec1_B, Vec vec2_B, Vec vec1_D, Vec vec2_D)
472: {
473:   PC_IS *pcis = (PC_IS *)(pc->data);

475:   if (!vec2_B) vec2_B = v;

477:   MatMult(pcis->A_BB, v, vec1_B);
478:   MatMult(pcis->A_IB, v, vec1_D);
479:   KSPSolve(pcis->ksp_D, vec1_D, vec2_D);
480:   KSPCheckSolve(pcis->ksp_D, pc, vec2_D);
481:   MatMult(pcis->A_BI, vec2_D, vec2_B);
482:   VecAXPY(vec1_B, -1.0, vec2_B);
483:   return 0;
484: }

486: /*
487:    PCISScatterArrayNToVecB - Scatters interface node values from a big array (of all local nodes, interior or interface,
488:    including ghosts) into an interface vector, when in SCATTER_FORWARD mode, or vice-versa, when in SCATTER_REVERSE
489:    mode.

491:    Input parameters:
492: .  pc - preconditioner context
493: .  array_N - [when in SCATTER_FORWARD mode] Array to be scattered into the vector
494: .  v_B - [when in SCATTER_REVERSE mode] Vector to be scattered into the array

496:    Output parameter:
497: .  array_N - [when in SCATTER_REVERSE mode] Array to receive the scattered vector
498: .  v_B - [when in SCATTER_FORWARD mode] Vector to receive the scattered array

500:    Note:
501:    The entries in the array that do not correspond to interface nodes remain unaltered.
502: */
503: PetscErrorCode PCISScatterArrayNToVecB(PetscScalar *array_N, Vec v_B, InsertMode imode, ScatterMode smode, PC pc)
504: {
505:   PetscInt        i;
506:   const PetscInt *idex;
507:   PetscScalar    *array_B;
508:   PC_IS          *pcis = (PC_IS *)(pc->data);

510:   VecGetArray(v_B, &array_B);
511:   ISGetIndices(pcis->is_B_local, &idex);

513:   if (smode == SCATTER_FORWARD) {
514:     if (imode == INSERT_VALUES) {
515:       for (i = 0; i < pcis->n_B; i++) array_B[i] = array_N[idex[i]];
516:     } else { /* ADD_VALUES */
517:       for (i = 0; i < pcis->n_B; i++) array_B[i] += array_N[idex[i]];
518:     }
519:   } else { /* SCATTER_REVERSE */
520:     if (imode == INSERT_VALUES) {
521:       for (i = 0; i < pcis->n_B; i++) array_N[idex[i]] = array_B[i];
522:     } else { /* ADD_VALUES */
523:       for (i = 0; i < pcis->n_B; i++) array_N[idex[i]] += array_B[i];
524:     }
525:   }
526:   ISRestoreIndices(pcis->is_B_local, &idex);
527:   VecRestoreArray(v_B, &array_B);
528:   return 0;
529: }

531: /*
532:    PCISApplyInvSchur - Solves the Neumann problem related to applying the inverse of the Schur complement.
533:    More precisely, solves the problem:
534:                                         [ A_II  A_IB ] [ . ]   [ 0 ]
535:                                         [            ] [   ] = [   ]
536:                                         [ A_BI  A_BB ] [ x ]   [ b ]

538:    Input parameters:
539: .  pc - preconditioner context
540: .  b - vector of local interface nodes (including ghosts)

542:    Output parameters:
543: .  x - vector of local interface nodes (including ghosts); returns the application of the inverse of the Schur
544:        complement to b
545: .  vec1_N - vector of local nodes (interior and interface, including ghosts); returns garbage (used as work space)
546: .  vec2_N - vector of local nodes (interior and interface, including ghosts); returns garbage (used as work space)

548: */
549: PetscErrorCode PCISApplyInvSchur(PC pc, Vec b, Vec x, Vec vec1_N, Vec vec2_N)
550: {
551:   PC_IS *pcis = (PC_IS *)(pc->data);

553:   /*
554:     Neumann solvers.
555:     Applying the inverse of the local Schur complement, i.e, solving a Neumann
556:     Problem with zero at the interior nodes of the RHS and extracting the interface
557:     part of the solution. inverse Schur complement is applied to b and the result
558:     is stored in x.
559:   */
560:   /* Setting the RHS vec1_N */
561:   VecSet(vec1_N, 0.0);
562:   VecScatterBegin(pcis->N_to_B, b, vec1_N, INSERT_VALUES, SCATTER_REVERSE);
563:   VecScatterEnd(pcis->N_to_B, b, vec1_N, INSERT_VALUES, SCATTER_REVERSE);
564:   /* Checking for consistency of the RHS */
565:   {
566:     PetscBool flg = PETSC_FALSE;
567:     PetscOptionsGetBool(NULL, NULL, "-pc_is_check_consistency", &flg, NULL);
568:     if (flg) {
569:       PetscScalar average;
570:       PetscViewer viewer;
571:       PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)pc), &viewer);

573:       VecSum(vec1_N, &average);
574:       average = average / ((PetscReal)pcis->n);
575:       PetscViewerASCIIPushSynchronized(viewer);
576:       if (pcis->pure_neumann) {
577:         PetscViewerASCIISynchronizedPrintf(viewer, "Subdomain %04d is floating. Average = % 1.14e\n", PetscGlobalRank, (double)PetscAbsScalar(average));
578:       } else {
579:         PetscViewerASCIISynchronizedPrintf(viewer, "Subdomain %04d is fixed.    Average = % 1.14e\n", PetscGlobalRank, (double)PetscAbsScalar(average));
580:       }
581:       PetscViewerFlush(viewer);
582:       PetscViewerASCIIPopSynchronized(viewer);
583:     }
584:   }
585:   /* Solving the system for vec2_N */
586:   KSPSolve(pcis->ksp_N, vec1_N, vec2_N);
587:   KSPCheckSolve(pcis->ksp_N, pc, vec2_N);
588:   /* Extracting the local interface vector out of the solution */
589:   VecScatterBegin(pcis->N_to_B, vec2_N, x, INSERT_VALUES, SCATTER_FORWARD);
590:   VecScatterEnd(pcis->N_to_B, vec2_N, x, INSERT_VALUES, SCATTER_FORWARD);
591:   return 0;
592: }