Actual source code: pcis.c

petsc-3.9.0 2018-04-07
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: {

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

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

 45:   PetscObjectReference((PetscObject)scaling_factors);
 46:   VecDestroy(&pcis->D);
 47:   pcis->D = scaling_factors;
 48:   return(0);
 49: }

 51: /*@
 52:  PCISSetSubdomainDiagonalScaling - Set diagonal scaling for PCIS.

 54:    Not collective

 56:    Input Parameters:
 57: +  pc - the preconditioning context
 58: -  scaling_factors - scaling factors for the subdomain

 60:    Level: intermediate

 62:    Notes:
 63:    Intended to use with jumping coefficients cases.

 65: .seealso: PCBDDC
 66: @*/
 67: PetscErrorCode PCISSetSubdomainDiagonalScaling(PC pc, Vec scaling_factors)
 68: {

 73:   PetscTryMethod(pc,"PCISSetSubdomainDiagonalScaling_C",(PC,Vec),(pc,scaling_factors));
 74:   return(0);
 75: }

 77: static PetscErrorCode PCISSetSubdomainScalingFactor_IS(PC pc, PetscScalar scal)
 78: {
 79:   PC_IS *pcis = (PC_IS*)pc->data;

 82:   pcis->scaling_factor = scal;
 83:   return(0);
 84: }

 86: /*@
 87:  PCISSetSubdomainScalingFactor - Set scaling factor for PCIS.

 89:    Not collective

 91:    Input Parameters:
 92: +  pc - the preconditioning context
 93: -  scal - scaling factor for the subdomain

 95:    Level: intermediate

 97:    Notes:
 98:    Intended to use with jumping coefficients cases.

100: .seealso: PCBDDC
101: @*/
102: PetscErrorCode PCISSetSubdomainScalingFactor(PC pc, PetscScalar scal)
103: {

108:   PetscTryMethod(pc,"PCISSetSubdomainScalingFactor_C",(PC,PetscScalar),(pc,scal));
109:   return(0);
110: }


113: /* -------------------------------------------------------------------------- */
114: /*
115:    PCISSetUp -
116: */
117: PetscErrorCode  PCISSetUp(PC pc, PetscBool computesolvers)
118: {
119:   PC_IS          *pcis  = (PC_IS*)(pc->data);
120:   Mat_IS         *matis;
121:   MatReuse       reuse;
123:   PetscBool      flg,issbaij;
124:   Vec            counter;

127:   PetscObjectTypeCompare((PetscObject)pc->pmat,MATIS,&flg);
128:   if (!flg) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONG,"Preconditioner type of Neumann Neumman requires matrix of type MATIS");
129:   matis = (Mat_IS*)pc->pmat->data;

131:   /* first time creation, get info on substructuring */
132:   if (!pc->setupcalled) {
133:     PetscInt    n_I;
134:     PetscInt    *idx_I_local,*idx_B_local,*idx_I_global,*idx_B_global;
135:     PetscBT     bt;
136:     PetscInt    i,j;

138:     /* get info on mapping */
139:     PetscObjectReference((PetscObject)pc->pmat->rmap->mapping);
140:     ISLocalToGlobalMappingDestroy(&pcis->mapping);
141:     pcis->mapping = pc->pmat->rmap->mapping;
142:     ISLocalToGlobalMappingGetSize(pcis->mapping,&pcis->n);
143:     ISLocalToGlobalMappingGetInfo(pcis->mapping,&(pcis->n_neigh),&(pcis->neigh),&(pcis->n_shared),&(pcis->shared));

145:     /* Identifying interior and interface nodes, in local numbering */
146:     PetscBTCreate(pcis->n,&bt);
147:     for (i=0;i<pcis->n_neigh;i++)
148:       for (j=0;j<pcis->n_shared[i];j++) {
149:           PetscBTSet(bt,pcis->shared[i][j]);
150:       }

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

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

171:     /* Creating the index sets */
172:     ISCreateGeneral(PETSC_COMM_SELF,pcis->n_B,idx_B_local,PETSC_COPY_VALUES, &pcis->is_B_local);
173:     ISCreateGeneral(PETSC_COMM_SELF,pcis->n_B,idx_B_global,PETSC_COPY_VALUES,&pcis->is_B_global);
174:     ISCreateGeneral(PETSC_COMM_SELF,n_I,idx_I_local,PETSC_COPY_VALUES, &pcis->is_I_local);
175:     ISCreateGeneral(PETSC_COMM_SELF,n_I,idx_I_global,PETSC_COPY_VALUES,&pcis->is_I_global);

177:     /* Freeing memory */
178:     PetscFree(idx_B_local);
179:     PetscFree(idx_I_local);
180:     PetscBTDestroy(&bt);

182:     /* Creating work vectors and arrays */
183:     VecDuplicate(matis->x,&pcis->vec1_N);
184:     VecDuplicate(pcis->vec1_N,&pcis->vec2_N);
185:     VecCreateSeq(PETSC_COMM_SELF,pcis->n-pcis->n_B,&pcis->vec1_D);
186:     VecDuplicate(pcis->vec1_D,&pcis->vec2_D);
187:     VecDuplicate(pcis->vec1_D,&pcis->vec3_D);
188:     VecDuplicate(pcis->vec1_D,&pcis->vec4_D);
189:     VecCreateSeq(PETSC_COMM_SELF,pcis->n_B,&pcis->vec1_B);
190:     VecDuplicate(pcis->vec1_B,&pcis->vec2_B);
191:     VecDuplicate(pcis->vec1_B,&pcis->vec3_B);
192:     MatCreateVecs(pc->pmat,&pcis->vec1_global,0);
193:     PetscMalloc1(pcis->n,&pcis->work_N);
194:     /* scaling vector */
195:     if (!pcis->D) { /* it can happen that the user passed in a scaling vector via PCISSetSubdomainDiagonalScaling */
196:       VecDuplicate(pcis->vec1_B,&pcis->D);
197:       VecSet(pcis->D,pcis->scaling_factor);
198:     }

200:     /* Creating the scatter contexts */
201:     VecScatterCreate(pcis->vec1_N,pcis->is_I_local,pcis->vec1_D,(IS)0,&pcis->N_to_D);
202:     VecScatterCreate(pcis->vec1_global,pcis->is_I_global,pcis->vec1_D,(IS)0,&pcis->global_to_D);
203:     VecScatterCreate(pcis->vec1_N,pcis->is_B_local,pcis->vec1_B,(IS)0,&pcis->N_to_B);
204:     VecScatterCreate(pcis->vec1_global,pcis->is_B_global,pcis->vec1_B,(IS)0,&pcis->global_to_B);

206:     /* map from boundary to local */
207:     ISLocalToGlobalMappingCreateIS(pcis->is_B_local,&pcis->BtoNmap);
208:   }

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

214:         [ A_II | A_IB ]
215:     A = [------+------]
216:         [ A_BI | A_BB ]
217:   */
218:   reuse = MAT_INITIAL_MATRIX;
219:   if (pcis->reusesubmatrices && pc->setupcalled) {
220:     if (pc->flag == SAME_NONZERO_PATTERN) {
221:       reuse = MAT_REUSE_MATRIX;
222:     } else {
223:       reuse = MAT_INITIAL_MATRIX;
224:     }
225:   }
226:   if (reuse == MAT_INITIAL_MATRIX) {
227:     MatDestroy(&pcis->A_II);
228:     MatDestroy(&pcis->A_IB);
229:     MatDestroy(&pcis->A_BI);
230:     MatDestroy(&pcis->A_BB);
231:   }

233:   MatCreateSubMatrix(matis->A,pcis->is_I_local,pcis->is_I_local,reuse,&pcis->A_II);
234:   MatCreateSubMatrix(matis->A,pcis->is_B_local,pcis->is_B_local,reuse,&pcis->A_BB);
235:   PetscObjectTypeCompare((PetscObject)matis->A,MATSEQSBAIJ,&issbaij);
236:   if (!issbaij) {
237:     MatCreateSubMatrix(matis->A,pcis->is_I_local,pcis->is_B_local,reuse,&pcis->A_IB);
238:     MatCreateSubMatrix(matis->A,pcis->is_B_local,pcis->is_I_local,reuse,&pcis->A_BI);
239:   } else {
240:     Mat newmat;
241:     MatConvert(matis->A,MATSEQBAIJ,MAT_INITIAL_MATRIX,&newmat);
242:     MatCreateSubMatrix(newmat,pcis->is_I_local,pcis->is_B_local,reuse,&pcis->A_IB);
243:     MatCreateSubMatrix(newmat,pcis->is_B_local,pcis->is_I_local,reuse,&pcis->A_BI);
244:     MatDestroy(&newmat);
245:   }

247:   /* Creating scaling vector D */
248:   PetscOptionsGetBool(((PetscObject)pc)->options,((PetscObject)pc)->prefix,"-pc_is_use_stiffness_scaling",&pcis->use_stiffness_scaling,NULL);
249:   if (pcis->use_stiffness_scaling) {
250:     PetscScalar *a;
251:     PetscInt    i,n;

253:     MatGetDiagonal(pcis->A_BB,pcis->D);
254:     VecGetLocalSize(pcis->D,&n);
255:     VecGetArray(pcis->D,&a);
256:     for (i=0;i<n;i++) if (PetscAbsScalar(a[i])<PETSC_SMALL) a[i] = 1.0;
257:     VecRestoreArray(pcis->D,&a);
258:   }
259:   MatCreateVecs(pc->pmat,&counter,0); /* temporary auxiliar vector */
260:   VecSet(counter,0.0);
261:   VecScatterBegin(pcis->global_to_B,pcis->D,counter,ADD_VALUES,SCATTER_REVERSE);
262:   VecScatterEnd(pcis->global_to_B,pcis->D,counter,ADD_VALUES,SCATTER_REVERSE);
263:   VecScatterBegin(pcis->global_to_B,counter,pcis->vec1_B,INSERT_VALUES,SCATTER_FORWARD);
264:   VecScatterEnd(pcis->global_to_B,counter,pcis->vec1_B,INSERT_VALUES,SCATTER_FORWARD);
265:   VecPointwiseDivide(pcis->D,pcis->D,pcis->vec1_B);
266:   VecDestroy(&counter);

268:   /* See historical note 01, at the bottom of this file. */

270:   /* Creating the KSP contexts for the local Dirichlet and Neumann problems */
271:   if (computesolvers) {
272:     PC pc_ctx;

274:     pcis->pure_neumann = matis->pure_neumann;
275:     /* Dirichlet */
276:     KSPCreate(PETSC_COMM_SELF,&pcis->ksp_D);
277:     KSPSetErrorIfNotConverged(pcis->ksp_D,pc->erroriffailure);
278:     PetscObjectIncrementTabLevel((PetscObject)pcis->ksp_D,(PetscObject)pc,1);
279:     KSPSetOperators(pcis->ksp_D,pcis->A_II,pcis->A_II);
280:     KSPSetOptionsPrefix(pcis->ksp_D,"is_localD_");
281:     KSPGetPC(pcis->ksp_D,&pc_ctx);
282:     PCSetType(pc_ctx,PCLU);
283:     KSPSetType(pcis->ksp_D,KSPPREONLY);
284:     KSPSetFromOptions(pcis->ksp_D);
285:     /* the vectors in the following line are dummy arguments, just telling the KSP the vector size. Values are not used */
286:     KSPSetUp(pcis->ksp_D);
287:     /* Neumann */
288:     KSPCreate(PETSC_COMM_SELF,&pcis->ksp_N);
289:     KSPSetErrorIfNotConverged(pcis->ksp_N,pc->erroriffailure);
290:     PetscObjectIncrementTabLevel((PetscObject)pcis->ksp_N,(PetscObject)pc,1);
291:     KSPSetOperators(pcis->ksp_N,matis->A,matis->A);
292:     KSPSetOptionsPrefix(pcis->ksp_N,"is_localN_");
293:     KSPGetPC(pcis->ksp_N,&pc_ctx);
294:     PCSetType(pc_ctx,PCLU);
295:     KSPSetType(pcis->ksp_N,KSPPREONLY);
296:     KSPSetFromOptions(pcis->ksp_N);
297:     {
298:       PetscBool damp_fixed                    = PETSC_FALSE,
299:                 remove_nullspace_fixed        = PETSC_FALSE,
300:                 set_damping_factor_floating   = PETSC_FALSE,
301:                 not_damp_floating             = PETSC_FALSE,
302:                 not_remove_nullspace_floating = PETSC_FALSE;
303:       PetscReal fixed_factor,
304:                 floating_factor;

306:       PetscOptionsGetReal(((PetscObject)pc_ctx)->options,((PetscObject)pc_ctx)->prefix,"-pc_is_damp_fixed",&fixed_factor,&damp_fixed);
307:       if (!damp_fixed) fixed_factor = 0.0;
308:       PetscOptionsGetBool(((PetscObject)pc_ctx)->options,((PetscObject)pc_ctx)->prefix,"-pc_is_damp_fixed",&damp_fixed,NULL);

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

312:       PetscOptionsGetReal(((PetscObject)pc_ctx)->options,((PetscObject)pc_ctx)->prefix,"-pc_is_set_damping_factor_floating",
313:                               &floating_factor,&set_damping_factor_floating);
314:       if (!set_damping_factor_floating) floating_factor = 0.0;
315:       PetscOptionsGetBool(((PetscObject)pc_ctx)->options,((PetscObject)pc_ctx)->prefix,"-pc_is_set_damping_factor_floating",&set_damping_factor_floating,NULL);
316:       if (!set_damping_factor_floating) floating_factor = 1.e-12;

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

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

322:       if (pcis->pure_neumann) {  /* floating subdomain */
323:         if (!(not_damp_floating)) {
324:           PCFactorSetShiftType(pc_ctx,MAT_SHIFT_NONZERO);
325:           PCFactorSetShiftAmount(pc_ctx,floating_factor);
326:         }
327:         if (!(not_remove_nullspace_floating)) {
328:           MatNullSpace nullsp;
329:           MatNullSpaceCreate(PETSC_COMM_SELF,PETSC_TRUE,0,NULL,&nullsp);
330:           MatSetNullSpace(matis->A,nullsp);
331:           MatNullSpaceDestroy(&nullsp);
332:         }
333:       } else {  /* fixed subdomain */
334:         if (damp_fixed) {
335:           PCFactorSetShiftType(pc_ctx,MAT_SHIFT_NONZERO);
336:           PCFactorSetShiftAmount(pc_ctx,floating_factor);
337:         }
338:         if (remove_nullspace_fixed) {
339:           MatNullSpace nullsp;
340:           MatNullSpaceCreate(PETSC_COMM_SELF,PETSC_TRUE,0,NULL,&nullsp);
341:           MatSetNullSpace(matis->A,nullsp);
342:           MatNullSpaceDestroy(&nullsp);
343:         }
344:       }
345:     }
346:     /* the vectors in the following line are dummy arguments, just telling the KSP the vector size. Values are not used */
347:     KSPSetUp(pcis->ksp_N);
348:   }

350:   return(0);
351: }

353: /* -------------------------------------------------------------------------- */
354: /*
355:    PCISDestroy -
356: */
357: PetscErrorCode  PCISDestroy(PC pc)
358: {
359:   PC_IS          *pcis = (PC_IS*)(pc->data);

363:   ISDestroy(&pcis->is_B_local);
364:   ISDestroy(&pcis->is_I_local);
365:   ISDestroy(&pcis->is_B_global);
366:   ISDestroy(&pcis->is_I_global);
367:   MatDestroy(&pcis->A_II);
368:   MatDestroy(&pcis->A_IB);
369:   MatDestroy(&pcis->A_BI);
370:   MatDestroy(&pcis->A_BB);
371:   VecDestroy(&pcis->D);
372:   KSPDestroy(&pcis->ksp_N);
373:   KSPDestroy(&pcis->ksp_D);
374:   VecDestroy(&pcis->vec1_N);
375:   VecDestroy(&pcis->vec2_N);
376:   VecDestroy(&pcis->vec1_D);
377:   VecDestroy(&pcis->vec2_D);
378:   VecDestroy(&pcis->vec3_D);
379:   VecDestroy(&pcis->vec4_D);
380:   VecDestroy(&pcis->vec1_B);
381:   VecDestroy(&pcis->vec2_B);
382:   VecDestroy(&pcis->vec3_B);
383:   VecDestroy(&pcis->vec1_global);
384:   VecScatterDestroy(&pcis->global_to_D);
385:   VecScatterDestroy(&pcis->N_to_B);
386:   VecScatterDestroy(&pcis->N_to_D);
387:   VecScatterDestroy(&pcis->global_to_B);
388:   PetscFree(pcis->work_N);
389:   if (pcis->n_neigh > -1) {
390:     ISLocalToGlobalMappingRestoreInfo(pcis->mapping,&(pcis->n_neigh),&(pcis->neigh),&(pcis->n_shared),&(pcis->shared));
391:   }
392:   ISLocalToGlobalMappingDestroy(&pcis->mapping);
393:   ISLocalToGlobalMappingDestroy(&pcis->BtoNmap);
394:   PetscObjectComposeFunction((PetscObject)pc,"PCISSetUseStiffnessScaling_C",NULL);
395:   PetscObjectComposeFunction((PetscObject)pc,"PCISSetSubdomainScalingFactor_C",NULL);
396:   PetscObjectComposeFunction((PetscObject)pc,"PCISSetSubdomainDiagonalScaling_C",NULL);
397:   return(0);
398: }

400: /* -------------------------------------------------------------------------- */
401: /*
402:    PCISCreate -
403: */
404: PetscErrorCode  PCISCreate(PC pc)
405: {
406:   PC_IS          *pcis = (PC_IS*)(pc->data);

410:   pcis->is_B_local       = 0;
411:   pcis->is_I_local       = 0;
412:   pcis->is_B_global      = 0;
413:   pcis->is_I_global      = 0;
414:   pcis->A_II             = 0;
415:   pcis->A_IB             = 0;
416:   pcis->A_BI             = 0;
417:   pcis->A_BB             = 0;
418:   pcis->D                = 0;
419:   pcis->ksp_N            = 0;
420:   pcis->ksp_D            = 0;
421:   pcis->vec1_N           = 0;
422:   pcis->vec2_N           = 0;
423:   pcis->vec1_D           = 0;
424:   pcis->vec2_D           = 0;
425:   pcis->vec3_D           = 0;
426:   pcis->vec1_B           = 0;
427:   pcis->vec2_B           = 0;
428:   pcis->vec3_B           = 0;
429:   pcis->vec1_global      = 0;
430:   pcis->work_N           = 0;
431:   pcis->global_to_D      = 0;
432:   pcis->N_to_B           = 0;
433:   pcis->N_to_D           = 0;
434:   pcis->global_to_B      = 0;
435:   pcis->mapping          = 0;
436:   pcis->BtoNmap          = 0;
437:   pcis->n_neigh          = -1;
438:   pcis->scaling_factor   = 1.0;
439:   pcis->reusesubmatrices = PETSC_TRUE;
440:   /* composing functions */
441:   PetscObjectComposeFunction((PetscObject)pc,"PCISSetUseStiffnessScaling_C",PCISSetUseStiffnessScaling_IS);
442:   PetscObjectComposeFunction((PetscObject)pc,"PCISSetSubdomainScalingFactor_C",PCISSetSubdomainScalingFactor_IS);
443:   PetscObjectComposeFunction((PetscObject)pc,"PCISSetSubdomainDiagonalScaling_C",PCISSetSubdomainDiagonalScaling_IS);
444:   return(0);
445: }

447: /* -------------------------------------------------------------------------- */
448: /*
449:    PCISApplySchur -

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

455:    Output parameters:
456: .  vec1_B - result of Schur complement applied to chunk
457: .  vec2_B - garbage (used as work space), or null (and v is used as workspace)
458: .  vec1_D - garbage (used as work space)
459: .  vec2_D - garbage (used as work space)

461: */
462: PetscErrorCode  PCISApplySchur(PC pc, Vec v, Vec vec1_B, Vec vec2_B, Vec vec1_D, Vec vec2_D)
463: {
465:   PC_IS          *pcis = (PC_IS*)(pc->data);

468:   if (!vec2_B) vec2_B = v;

470:   MatMult(pcis->A_BB,v,vec1_B);
471:   MatMult(pcis->A_IB,v,vec1_D);
472:   KSPSolve(pcis->ksp_D,vec1_D,vec2_D);
473:   MatMult(pcis->A_BI,vec2_D,vec2_B);
474:   VecAXPY(vec1_B,-1.0,vec2_B);
475:   return(0);
476: }

478: /* -------------------------------------------------------------------------- */
479: /*
480:    PCISScatterArrayNToVecB - Scatters interface node values from a big array (of all local nodes, interior or interface,
481:    including ghosts) into an interface vector, when in SCATTER_FORWARD mode, or vice-versa, when in SCATTER_REVERSE
482:    mode.

484:    Input parameters:
485: .  pc - preconditioner context
486: .  array_N - [when in SCATTER_FORWARD mode] Array to be scattered into the vector
487: .  v_B - [when in SCATTER_REVERSE mode] Vector to be scattered into the array

489:    Output parameter:
490: .  array_N - [when in SCATTER_REVERSE mode] Array to receive the scattered vector
491: .  v_B - [when in SCATTER_FORWARD mode] Vector to receive the scattered array

493:    Notes:
494:    The entries in the array that do not correspond to interface nodes remain unaltered.
495: */
496: PetscErrorCode  PCISScatterArrayNToVecB(PetscScalar *array_N, Vec v_B, InsertMode imode, ScatterMode smode, PC pc)
497: {
498:   PetscInt       i;
499:   const PetscInt *idex;
501:   PetscScalar    *array_B;
502:   PC_IS          *pcis = (PC_IS*)(pc->data);

505:   VecGetArray(v_B,&array_B);
506:   ISGetIndices(pcis->is_B_local,&idex);

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

526: /* -------------------------------------------------------------------------- */
527: /*
528:    PCISApplyInvSchur - Solves the Neumann problem related to applying the inverse of the Schur complement.
529:    More precisely, solves the problem:
530:                                         [ A_II  A_IB ] [ . ]   [ 0 ]
531:                                         [            ] [   ] = [   ]
532:                                         [ A_BI  A_BB ] [ x ]   [ b ]

534:    Input parameters:
535: .  pc - preconditioner context
536: .  b - vector of local interface nodes (including ghosts)

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

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

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

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