Actual source code: pcis.c

petsc-3.14.6 2021-03-30
Report Typos and Errors

  2: #include <../src/ksp/pc/impls/is/pcis.h>

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

  9:   pcis->use_stiffness_scaling = use;
 10:   return(0);
 11: }

 13: /*@
 14:  PCISSetUseStiffnessScaling - Tells PCIS to construct partition of unity using
 15:                               local matrices' diagonal.

 17:    Not collective

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

 23:    Level: intermediate

 25:    Notes:

 27: .seealso: PCBDDC
 28: @*/
 29: PetscErrorCode PCISSetUseStiffnessScaling(PC pc, PetscBool use)
 30: {

 36:   PetscTryMethod(pc,"PCISSetUseStiffnessScaling_C",(PC,PetscBool),(pc,use));
 37:   return(0);
 38: }

 40: static PetscErrorCode PCISSetSubdomainDiagonalScaling_IS(PC pc, Vec scaling_factors)
 41: {
 43:   PC_IS          *pcis = (PC_IS*)pc->data;

 46:   PetscObjectReference((PetscObject)scaling_factors);
 47:   VecDestroy(&pcis->D);
 48:   pcis->D = scaling_factors;
 49:   if (pc->setupcalled) {
 50:     PetscInt sn;

 52:     VecGetSize(pcis->D,&sn);
 53:     if (sn == pcis->n) {
 54:       VecScatterBegin(pcis->N_to_B,pcis->D,pcis->vec1_B,INSERT_VALUES,SCATTER_FORWARD);
 55:       VecScatterEnd(pcis->N_to_B,pcis->D,pcis->vec1_B,INSERT_VALUES,SCATTER_FORWARD);
 56:       VecDestroy(&pcis->D);
 57:       VecDuplicate(pcis->vec1_B,&pcis->D);
 58:       VecCopy(pcis->vec1_B,pcis->D);
 59:     } else if (sn != pcis->n_B) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Invalid size for scaling vector. Expected %D (or full %D), found %D",pcis->n_B,pcis->n,sn);
 60:   }
 61:   return(0);
 62: }

 64: /*@
 65:  PCISSetSubdomainDiagonalScaling - Set diagonal scaling for PCIS.

 67:    Not collective

 69:    Input Parameters:
 70: +  pc - the preconditioning context
 71: -  scaling_factors - scaling factors for the subdomain

 73:    Level: intermediate

 75:    Notes:
 76:    Intended to use with jumping coefficients cases.

 78: .seealso: PCBDDC
 79: @*/
 80: PetscErrorCode PCISSetSubdomainDiagonalScaling(PC pc, Vec scaling_factors)
 81: {

 87:   PetscTryMethod(pc,"PCISSetSubdomainDiagonalScaling_C",(PC,Vec),(pc,scaling_factors));
 88:   return(0);
 89: }

 91: static PetscErrorCode PCISSetSubdomainScalingFactor_IS(PC pc, PetscScalar scal)
 92: {
 93:   PC_IS *pcis = (PC_IS*)pc->data;

 96:   pcis->scaling_factor = scal;
 97:   if (pcis->D) {

100:     VecSet(pcis->D,pcis->scaling_factor);
101:   }
102:   return(0);
103: }

105: /*@
106:  PCISSetSubdomainScalingFactor - Set scaling factor for PCIS.

108:    Not collective

110:    Input Parameters:
111: +  pc - the preconditioning context
112: -  scal - scaling factor for the subdomain

114:    Level: intermediate

116:    Notes:
117:    Intended to use with jumping coefficients cases.

119: .seealso: PCBDDC
120: @*/
121: PetscErrorCode PCISSetSubdomainScalingFactor(PC pc, PetscScalar scal)
122: {

127:   PetscTryMethod(pc,"PCISSetSubdomainScalingFactor_C",(PC,PetscScalar),(pc,scal));
128:   return(0);
129: }


132: /* -------------------------------------------------------------------------- */
133: /*
134:    PCISSetUp -
135: */
136: PetscErrorCode  PCISSetUp(PC pc, PetscBool computematrices, PetscBool computesolvers)
137: {
138:   PC_IS          *pcis  = (PC_IS*)(pc->data);
139:   Mat_IS         *matis;
140:   MatReuse       reuse;
142:   PetscBool      flg,issbaij;

145:   PetscObjectTypeCompare((PetscObject)pc->pmat,MATIS,&flg);
146:   if (!flg) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONG,"Requires preconditioning matrix of type MATIS");
147:   matis = (Mat_IS*)pc->pmat->data;
148:   if (pc->useAmat) {
149:     PetscObjectTypeCompare((PetscObject)pc->mat,MATIS,&flg);
150:     if (!flg) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONG,"Requires linear system matrix of type MATIS");
151:   }

153:   /* first time creation, get info on substructuring */
154:   if (!pc->setupcalled) {
155:     PetscInt    n_I;
156:     PetscInt    *idx_I_local,*idx_B_local,*idx_I_global,*idx_B_global;
157:     PetscBT     bt;
158:     PetscInt    i,j;

160:     /* get info on mapping */
161:     PetscObjectReference((PetscObject)pc->pmat->rmap->mapping);
162:     ISLocalToGlobalMappingDestroy(&pcis->mapping);
163:     pcis->mapping = pc->pmat->rmap->mapping;
164:     ISLocalToGlobalMappingGetSize(pcis->mapping,&pcis->n);
165:     ISLocalToGlobalMappingGetInfo(pcis->mapping,&(pcis->n_neigh),&(pcis->neigh),&(pcis->n_shared),&(pcis->shared));

167:     /* Identifying interior and interface nodes, in local numbering */
168:     PetscBTCreate(pcis->n,&bt);
169:     for (i=0;i<pcis->n_neigh;i++)
170:       for (j=0;j<pcis->n_shared[i];j++) {
171:         PetscBTSet(bt,pcis->shared[i][j]);
172:       }

174:     /* Creating local and global index sets for interior and inteface nodes. */
175:     PetscMalloc1(pcis->n,&idx_I_local);
176:     PetscMalloc1(pcis->n,&idx_B_local);
177:     for (i=0, pcis->n_B=0, n_I=0; i<pcis->n; i++) {
178:       if (!PetscBTLookup(bt,i)) {
179:         idx_I_local[n_I] = i;
180:         n_I++;
181:       } else {
182:         idx_B_local[pcis->n_B] = i;
183:         pcis->n_B++;
184:       }
185:     }

187:     /* Getting the global numbering */
188:     idx_B_global = idx_I_local + n_I; /* Just avoiding allocating extra memory, since we have vacant space */
189:     idx_I_global = idx_B_local + pcis->n_B;
190:     ISLocalToGlobalMappingApply(pcis->mapping,pcis->n_B,idx_B_local,idx_B_global);
191:     ISLocalToGlobalMappingApply(pcis->mapping,n_I,idx_I_local,idx_I_global);

193:     /* Creating the index sets */
194:     ISCreateGeneral(PETSC_COMM_SELF,pcis->n_B,idx_B_local,PETSC_COPY_VALUES, &pcis->is_B_local);
195:     ISCreateGeneral(PetscObjectComm((PetscObject)pc),pcis->n_B,idx_B_global,PETSC_COPY_VALUES,&pcis->is_B_global);
196:     ISCreateGeneral(PETSC_COMM_SELF,n_I,idx_I_local,PETSC_COPY_VALUES, &pcis->is_I_local);
197:     ISCreateGeneral(PetscObjectComm((PetscObject)pc),n_I,idx_I_global,PETSC_COPY_VALUES,&pcis->is_I_global);

199:     /* Freeing memory */
200:     PetscFree(idx_B_local);
201:     PetscFree(idx_I_local);
202:     PetscBTDestroy(&bt);

204:     /* Creating work vectors and arrays */
205:     VecDuplicate(matis->x,&pcis->vec1_N);
206:     VecDuplicate(pcis->vec1_N,&pcis->vec2_N);
207:     VecCreate(PETSC_COMM_SELF,&pcis->vec1_D);
208:     VecSetSizes(pcis->vec1_D,pcis->n-pcis->n_B,PETSC_DECIDE);
209:     VecSetType(pcis->vec1_D,((PetscObject)pcis->vec1_N)->type_name);
210:     VecDuplicate(pcis->vec1_D,&pcis->vec2_D);
211:     VecDuplicate(pcis->vec1_D,&pcis->vec3_D);
212:     VecDuplicate(pcis->vec1_D,&pcis->vec4_D);
213:     VecCreate(PETSC_COMM_SELF,&pcis->vec1_B);
214:     VecSetSizes(pcis->vec1_B,pcis->n_B,PETSC_DECIDE);
215:     VecSetType(pcis->vec1_B,((PetscObject)pcis->vec1_N)->type_name);
216:     VecDuplicate(pcis->vec1_B,&pcis->vec2_B);
217:     VecDuplicate(pcis->vec1_B,&pcis->vec3_B);
218:     MatCreateVecs(pc->pmat,&pcis->vec1_global,NULL);
219:     PetscMalloc1(pcis->n,&pcis->work_N);
220:     /* scaling vector */
221:     if (!pcis->D) { /* it can happen that the user passed in a scaling vector via PCISSetSubdomainDiagonalScaling */
222:       VecDuplicate(pcis->vec1_B,&pcis->D);
223:       VecSet(pcis->D,pcis->scaling_factor);
224:     }

226:     /* Creating the scatter contexts */
227:     VecScatterCreate(pcis->vec1_N,pcis->is_I_local,pcis->vec1_D,(IS)0,&pcis->N_to_D);
228:     VecScatterCreate(pcis->vec1_global,pcis->is_I_global,pcis->vec1_D,(IS)0,&pcis->global_to_D);
229:     VecScatterCreate(pcis->vec1_N,pcis->is_B_local,pcis->vec1_B,(IS)0,&pcis->N_to_B);
230:     VecScatterCreate(pcis->vec1_global,pcis->is_B_global,pcis->vec1_B,(IS)0,&pcis->global_to_B);

232:     /* map from boundary to local */
233:     ISLocalToGlobalMappingCreateIS(pcis->is_B_local,&pcis->BtoNmap);
234:   }

236:   {
237:     PetscInt sn;

239:     VecGetSize(pcis->D,&sn);
240:     if (sn == pcis->n) {
241:       VecScatterBegin(pcis->N_to_B,pcis->D,pcis->vec1_B,INSERT_VALUES,SCATTER_FORWARD);
242:       VecScatterEnd(pcis->N_to_B,pcis->D,pcis->vec1_B,INSERT_VALUES,SCATTER_FORWARD);
243:       VecDestroy(&pcis->D);
244:       VecDuplicate(pcis->vec1_B,&pcis->D);
245:       VecCopy(pcis->vec1_B,pcis->D);
246:     } else if (sn != pcis->n_B) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Invalid size for scaling vector. Expected %D (or full %D), found %D",pcis->n_B,pcis->n,sn);
247:   }

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

253:         [ A_II | A_IB ]
254:     A = [------+------]
255:         [ A_BI | A_BB ]
256:   */
257:   if (computematrices) {
258:     PetscBool amat = (PetscBool)(pc->mat != pc->pmat && pc->useAmat);
259:     PetscInt  bs,ibs;

261:     reuse = MAT_INITIAL_MATRIX;
262:     if (pcis->reusesubmatrices && pc->setupcalled) {
263:       if (pc->flag == SAME_NONZERO_PATTERN) {
264:         reuse = MAT_REUSE_MATRIX;
265:       } else {
266:         reuse = MAT_INITIAL_MATRIX;
267:       }
268:     }
269:     if (reuse == MAT_INITIAL_MATRIX) {
270:       MatDestroy(&pcis->A_II);
271:       MatDestroy(&pcis->pA_II);
272:       MatDestroy(&pcis->A_IB);
273:       MatDestroy(&pcis->A_BI);
274:       MatDestroy(&pcis->A_BB);
275:     }

277:     ISLocalToGlobalMappingGetBlockSize(pcis->mapping,&ibs);
278:     MatGetBlockSize(matis->A,&bs);
279:     MatCreateSubMatrix(matis->A,pcis->is_I_local,pcis->is_I_local,reuse,&pcis->pA_II);
280:     if (amat) {
281:       Mat_IS *amatis = (Mat_IS*)pc->mat->data;
282:       MatCreateSubMatrix(amatis->A,pcis->is_I_local,pcis->is_I_local,reuse,&pcis->A_II);
283:     } else {
284:       PetscObjectReference((PetscObject)pcis->pA_II);
285:       MatDestroy(&pcis->A_II);
286:       pcis->A_II = pcis->pA_II;
287:     }
288:     MatSetBlockSize(pcis->A_II,bs == ibs ? bs : 1);
289:     MatSetBlockSize(pcis->pA_II,bs == ibs ? bs : 1);
290:     MatCreateSubMatrix(matis->A,pcis->is_B_local,pcis->is_B_local,reuse,&pcis->A_BB);
291:     PetscObjectTypeCompare((PetscObject)matis->A,MATSEQSBAIJ,&issbaij);
292:     if (!issbaij) {
293:       MatCreateSubMatrix(matis->A,pcis->is_I_local,pcis->is_B_local,reuse,&pcis->A_IB);
294:       MatCreateSubMatrix(matis->A,pcis->is_B_local,pcis->is_I_local,reuse,&pcis->A_BI);
295:     } else {
296:       Mat newmat;

298:       MatConvert(matis->A,MATSEQBAIJ,MAT_INITIAL_MATRIX,&newmat);
299:       MatCreateSubMatrix(newmat,pcis->is_I_local,pcis->is_B_local,reuse,&pcis->A_IB);
300:       MatCreateSubMatrix(newmat,pcis->is_B_local,pcis->is_I_local,reuse,&pcis->A_BI);
301:       MatDestroy(&newmat);
302:     }
303:     MatSetBlockSize(pcis->A_BB,bs == ibs ? bs : 1);
304:   }

306:   /* Creating scaling vector D */
307:   PetscOptionsGetBool(((PetscObject)pc)->options,((PetscObject)pc)->prefix,"-pc_is_use_stiffness_scaling",&pcis->use_stiffness_scaling,NULL);
308:   if (pcis->use_stiffness_scaling) {
309:     PetscScalar *a;
310:     PetscInt    i,n;

312:     if (pcis->A_BB) {
313:       MatGetDiagonal(pcis->A_BB,pcis->D);
314:     } else {
315:       MatGetDiagonal(matis->A,pcis->vec1_N);
316:       VecScatterBegin(pcis->N_to_B,pcis->vec1_N,pcis->D,INSERT_VALUES,SCATTER_FORWARD);
317:       VecScatterEnd(pcis->N_to_B,pcis->vec1_N,pcis->D,INSERT_VALUES,SCATTER_FORWARD);
318:     }
319:     VecAbs(pcis->D);
320:     VecGetLocalSize(pcis->D,&n);
321:     VecGetArray(pcis->D,&a);
322:     for (i=0;i<n;i++) if (PetscAbsScalar(a[i])<PETSC_SMALL) a[i] = 1.0;
323:     VecRestoreArray(pcis->D,&a);
324:   }
325:   VecSet(pcis->vec1_global,0.0);
326:   VecScatterBegin(pcis->global_to_B,pcis->D,pcis->vec1_global,ADD_VALUES,SCATTER_REVERSE);
327:   VecScatterEnd(pcis->global_to_B,pcis->D,pcis->vec1_global,ADD_VALUES,SCATTER_REVERSE);
328:   VecScatterBegin(pcis->global_to_B,pcis->vec1_global,pcis->vec1_B,INSERT_VALUES,SCATTER_FORWARD);
329:   VecScatterEnd(pcis->global_to_B,pcis->vec1_global,pcis->vec1_B,INSERT_VALUES,SCATTER_FORWARD);
330:   VecPointwiseDivide(pcis->D,pcis->D,pcis->vec1_B);
331:   /* See historical note 01, at the bottom of this file. */

333:   /* Creating the KSP contexts for the local Dirichlet and Neumann problems */
334:   if (computesolvers) {
335:     PC pc_ctx;

337:     pcis->pure_neumann = matis->pure_neumann;
338:     /* Dirichlet */
339:     KSPCreate(PETSC_COMM_SELF,&pcis->ksp_D);
340:     KSPSetErrorIfNotConverged(pcis->ksp_D,pc->erroriffailure);
341:     PetscObjectIncrementTabLevel((PetscObject)pcis->ksp_D,(PetscObject)pc,1);
342:     KSPSetOperators(pcis->ksp_D,pcis->A_II,pcis->A_II);
343:     KSPSetOptionsPrefix(pcis->ksp_D,"is_localD_");
344:     KSPGetPC(pcis->ksp_D,&pc_ctx);
345:     PCSetType(pc_ctx,PCLU);
346:     KSPSetType(pcis->ksp_D,KSPPREONLY);
347:     KSPSetFromOptions(pcis->ksp_D);
348:     /* the vectors in the following line are dummy arguments, just telling the KSP the vector size. Values are not used */
349:     KSPSetUp(pcis->ksp_D);
350:     /* Neumann */
351:     KSPCreate(PETSC_COMM_SELF,&pcis->ksp_N);
352:     KSPSetErrorIfNotConverged(pcis->ksp_N,pc->erroriffailure);
353:     PetscObjectIncrementTabLevel((PetscObject)pcis->ksp_N,(PetscObject)pc,1);
354:     KSPSetOperators(pcis->ksp_N,matis->A,matis->A);
355:     KSPSetOptionsPrefix(pcis->ksp_N,"is_localN_");
356:     KSPGetPC(pcis->ksp_N,&pc_ctx);
357:     PCSetType(pc_ctx,PCLU);
358:     KSPSetType(pcis->ksp_N,KSPPREONLY);
359:     KSPSetFromOptions(pcis->ksp_N);
360:     {
361:       PetscBool damp_fixed                    = PETSC_FALSE,
362:                 remove_nullspace_fixed        = PETSC_FALSE,
363:                 set_damping_factor_floating   = PETSC_FALSE,
364:                 not_damp_floating             = PETSC_FALSE,
365:                 not_remove_nullspace_floating = PETSC_FALSE;
366:       PetscReal fixed_factor,
367:                 floating_factor;

369:       PetscOptionsGetReal(((PetscObject)pc_ctx)->options,((PetscObject)pc_ctx)->prefix,"-pc_is_damp_fixed",&fixed_factor,&damp_fixed);
370:       if (!damp_fixed) fixed_factor = 0.0;
371:       PetscOptionsGetBool(((PetscObject)pc_ctx)->options,((PetscObject)pc_ctx)->prefix,"-pc_is_damp_fixed",&damp_fixed,NULL);

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

375:       PetscOptionsGetReal(((PetscObject)pc_ctx)->options,((PetscObject)pc_ctx)->prefix,"-pc_is_set_damping_factor_floating",
376:                               &floating_factor,&set_damping_factor_floating);
377:       if (!set_damping_factor_floating) floating_factor = 0.0;
378:       PetscOptionsGetBool(((PetscObject)pc_ctx)->options,((PetscObject)pc_ctx)->prefix,"-pc_is_set_damping_factor_floating",&set_damping_factor_floating,NULL);
379:       if (!set_damping_factor_floating) floating_factor = 1.e-12;

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

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

385:       if (pcis->pure_neumann) {  /* floating subdomain */
386:         if (!(not_damp_floating)) {
387:           PCFactorSetShiftType(pc_ctx,MAT_SHIFT_NONZERO);
388:           PCFactorSetShiftAmount(pc_ctx,floating_factor);
389:         }
390:         if (!(not_remove_nullspace_floating)) {
391:           MatNullSpace nullsp;
392:           MatNullSpaceCreate(PETSC_COMM_SELF,PETSC_TRUE,0,NULL,&nullsp);
393:           MatSetNullSpace(matis->A,nullsp);
394:           MatNullSpaceDestroy(&nullsp);
395:         }
396:       } else {  /* fixed subdomain */
397:         if (damp_fixed) {
398:           PCFactorSetShiftType(pc_ctx,MAT_SHIFT_NONZERO);
399:           PCFactorSetShiftAmount(pc_ctx,floating_factor);
400:         }
401:         if (remove_nullspace_fixed) {
402:           MatNullSpace nullsp;
403:           MatNullSpaceCreate(PETSC_COMM_SELF,PETSC_TRUE,0,NULL,&nullsp);
404:           MatSetNullSpace(matis->A,nullsp);
405:           MatNullSpaceDestroy(&nullsp);
406:         }
407:       }
408:     }
409:     /* the vectors in the following line are dummy arguments, just telling the KSP the vector size. Values are not used */
410:     KSPSetUp(pcis->ksp_N);
411:   }
412:   return(0);
413: }

415: /* -------------------------------------------------------------------------- */
416: /*
417:    PCISDestroy -
418: */
419: PetscErrorCode  PCISDestroy(PC pc)
420: {
421:   PC_IS          *pcis = (PC_IS*)(pc->data);

425:   ISDestroy(&pcis->is_B_local);
426:   ISDestroy(&pcis->is_I_local);
427:   ISDestroy(&pcis->is_B_global);
428:   ISDestroy(&pcis->is_I_global);
429:   MatDestroy(&pcis->A_II);
430:   MatDestroy(&pcis->pA_II);
431:   MatDestroy(&pcis->A_IB);
432:   MatDestroy(&pcis->A_BI);
433:   MatDestroy(&pcis->A_BB);
434:   VecDestroy(&pcis->D);
435:   KSPDestroy(&pcis->ksp_N);
436:   KSPDestroy(&pcis->ksp_D);
437:   VecDestroy(&pcis->vec1_N);
438:   VecDestroy(&pcis->vec2_N);
439:   VecDestroy(&pcis->vec1_D);
440:   VecDestroy(&pcis->vec2_D);
441:   VecDestroy(&pcis->vec3_D);
442:   VecDestroy(&pcis->vec4_D);
443:   VecDestroy(&pcis->vec1_B);
444:   VecDestroy(&pcis->vec2_B);
445:   VecDestroy(&pcis->vec3_B);
446:   VecDestroy(&pcis->vec1_global);
447:   VecScatterDestroy(&pcis->global_to_D);
448:   VecScatterDestroy(&pcis->N_to_B);
449:   VecScatterDestroy(&pcis->N_to_D);
450:   VecScatterDestroy(&pcis->global_to_B);
451:   PetscFree(pcis->work_N);
452:   if (pcis->n_neigh > -1) {
453:     ISLocalToGlobalMappingRestoreInfo(pcis->mapping,&(pcis->n_neigh),&(pcis->neigh),&(pcis->n_shared),&(pcis->shared));
454:   }
455:   ISLocalToGlobalMappingDestroy(&pcis->mapping);
456:   ISLocalToGlobalMappingDestroy(&pcis->BtoNmap);
457:   PetscObjectComposeFunction((PetscObject)pc,"PCISSetUseStiffnessScaling_C",NULL);
458:   PetscObjectComposeFunction((PetscObject)pc,"PCISSetSubdomainScalingFactor_C",NULL);
459:   PetscObjectComposeFunction((PetscObject)pc,"PCISSetSubdomainDiagonalScaling_C",NULL);
460:   return(0);
461: }

463: /* -------------------------------------------------------------------------- */
464: /*
465:    PCISCreate -
466: */
467: PetscErrorCode  PCISCreate(PC pc)
468: {
469:   PC_IS          *pcis = (PC_IS*)(pc->data);

473:   pcis->n_neigh          = -1;
474:   pcis->scaling_factor   = 1.0;
475:   pcis->reusesubmatrices = PETSC_TRUE;
476:   /* composing functions */
477:   PetscObjectComposeFunction((PetscObject)pc,"PCISSetUseStiffnessScaling_C",PCISSetUseStiffnessScaling_IS);
478:   PetscObjectComposeFunction((PetscObject)pc,"PCISSetSubdomainScalingFactor_C",PCISSetSubdomainScalingFactor_IS);
479:   PetscObjectComposeFunction((PetscObject)pc,"PCISSetSubdomainDiagonalScaling_C",PCISSetSubdomainDiagonalScaling_IS);
480:   return(0);
481: }

483: /* -------------------------------------------------------------------------- */
484: /*
485:    PCISApplySchur -

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

491:    Output parameters:
492: .  vec1_B - result of Schur complement applied to chunk
493: .  vec2_B - garbage (used as work space), or null (and v is used as workspace)
494: .  vec1_D - garbage (used as work space)
495: .  vec2_D - garbage (used as work space)

497: */
498: PetscErrorCode  PCISApplySchur(PC pc, Vec v, Vec vec1_B, Vec vec2_B, Vec vec1_D, Vec vec2_D)
499: {
501:   PC_IS          *pcis = (PC_IS*)(pc->data);

504:   if (!vec2_B) vec2_B = v;

506:   MatMult(pcis->A_BB,v,vec1_B);
507:   MatMult(pcis->A_IB,v,vec1_D);
508:   KSPSolve(pcis->ksp_D,vec1_D,vec2_D);
509:   KSPCheckSolve(pcis->ksp_D,pc,vec2_D);
510:   MatMult(pcis->A_BI,vec2_D,vec2_B);
511:   VecAXPY(vec1_B,-1.0,vec2_B);
512:   return(0);
513: }

515: /* -------------------------------------------------------------------------- */
516: /*
517:    PCISScatterArrayNToVecB - Scatters interface node values from a big array (of all local nodes, interior or interface,
518:    including ghosts) into an interface vector, when in SCATTER_FORWARD mode, or vice-versa, when in SCATTER_REVERSE
519:    mode.

521:    Input parameters:
522: .  pc - preconditioner context
523: .  array_N - [when in SCATTER_FORWARD mode] Array to be scattered into the vector
524: .  v_B - [when in SCATTER_REVERSE mode] Vector to be scattered into the array

526:    Output parameter:
527: .  array_N - [when in SCATTER_REVERSE mode] Array to receive the scattered vector
528: .  v_B - [when in SCATTER_FORWARD mode] Vector to receive the scattered array

530:    Notes:
531:    The entries in the array that do not correspond to interface nodes remain unaltered.
532: */
533: PetscErrorCode  PCISScatterArrayNToVecB(PetscScalar *array_N, Vec v_B, InsertMode imode, ScatterMode smode, PC pc)
534: {
535:   PetscInt       i;
536:   const PetscInt *idex;
538:   PetscScalar    *array_B;
539:   PC_IS          *pcis = (PC_IS*)(pc->data);

542:   VecGetArray(v_B,&array_B);
543:   ISGetIndices(pcis->is_B_local,&idex);

545:   if (smode == SCATTER_FORWARD) {
546:     if (imode == INSERT_VALUES) {
547:       for (i=0; i<pcis->n_B; i++) array_B[i] = array_N[idex[i]];
548:     } else {  /* ADD_VALUES */
549:       for (i=0; i<pcis->n_B; i++) array_B[i] += array_N[idex[i]];
550:     }
551:   } else {  /* SCATTER_REVERSE */
552:     if (imode == INSERT_VALUES) {
553:       for (i=0; i<pcis->n_B; i++) array_N[idex[i]] = array_B[i];
554:     } else {  /* ADD_VALUES */
555:       for (i=0; i<pcis->n_B; i++) array_N[idex[i]] += array_B[i];
556:     }
557:   }
558:   ISRestoreIndices(pcis->is_B_local,&idex);
559:   VecRestoreArray(v_B,&array_B);
560:   return(0);
561: }

563: /* -------------------------------------------------------------------------- */
564: /*
565:    PCISApplyInvSchur - Solves the Neumann problem related to applying the inverse of the Schur complement.
566:    More precisely, solves the problem:
567:                                         [ A_II  A_IB ] [ . ]   [ 0 ]
568:                                         [            ] [   ] = [   ]
569:                                         [ A_BI  A_BB ] [ x ]   [ b ]

571:    Input parameters:
572: .  pc - preconditioner context
573: .  b - vector of local interface nodes (including ghosts)

575:    Output parameters:
576: .  x - vector of local interface nodes (including ghosts); returns the application of the inverse of the Schur
577:        complement to b
578: .  vec1_N - vector of local nodes (interior and interface, including ghosts); returns garbage (used as work space)
579: .  vec2_N - vector of local nodes (interior and interface, including ghosts); returns garbage (used as work space)

581: */
582: PetscErrorCode  PCISApplyInvSchur(PC pc, Vec b, Vec x, Vec vec1_N, Vec vec2_N)
583: {
585:   PC_IS          *pcis = (PC_IS*)(pc->data);

588:   /*
589:     Neumann solvers.
590:     Applying the inverse of the local Schur complement, i.e, solving a Neumann
591:     Problem with zero at the interior nodes of the RHS and extracting the interface
592:     part of the solution. inverse Schur complement is applied to b and the result
593:     is stored in x.
594:   */
595:   /* Setting the RHS vec1_N */
596:   VecSet(vec1_N,0.0);
597:   VecScatterBegin(pcis->N_to_B,b,vec1_N,INSERT_VALUES,SCATTER_REVERSE);
598:   VecScatterEnd  (pcis->N_to_B,b,vec1_N,INSERT_VALUES,SCATTER_REVERSE);
599:   /* Checking for consistency of the RHS */
600:   {
601:     PetscBool flg = PETSC_FALSE;
602:     PetscOptionsGetBool(NULL,NULL,"-pc_is_check_consistency",&flg,NULL);
603:     if (flg) {
604:       PetscScalar average;
605:       PetscViewer viewer;
606:       PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)pc),&viewer);

608:       VecSum(vec1_N,&average);
609:       average = average / ((PetscReal)pcis->n);
610:       PetscViewerASCIIPushSynchronized(viewer);
611:       if (pcis->pure_neumann) {
612:         PetscViewerASCIISynchronizedPrintf(viewer,"Subdomain %04d is floating. Average = % 1.14e\n",PetscGlobalRank,PetscAbsScalar(average));
613:       } else {
614:         PetscViewerASCIISynchronizedPrintf(viewer,"Subdomain %04d is fixed.    Average = % 1.14e\n",PetscGlobalRank,PetscAbsScalar(average));
615:       }
616:       PetscViewerFlush(viewer);
617:       PetscViewerASCIIPopSynchronized(viewer);
618:     }
619:   }
620:   /* Solving the system for vec2_N */
621:   KSPSolve(pcis->ksp_N,vec1_N,vec2_N);
622:   KSPCheckSolve(pcis->ksp_N,pc,vec2_N);
623:   /* Extracting the local interface vector out of the solution */
624:   VecScatterBegin(pcis->N_to_B,vec2_N,x,INSERT_VALUES,SCATTER_FORWARD);
625:   VecScatterEnd  (pcis->N_to_B,vec2_N,x,INSERT_VALUES,SCATTER_FORWARD);
626:   return(0);
627: }