Actual source code: pcis.c
petsc-3.8.4 2018-03-24
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: MatGetDiagonal(pcis->A_BB,pcis->D);
251: }
252: MatCreateVecs(pc->pmat,&counter,0); /* temporary auxiliar vector */
253: VecSet(counter,0.0);
254: VecScatterBegin(pcis->global_to_B,pcis->D,counter,ADD_VALUES,SCATTER_REVERSE);
255: VecScatterEnd(pcis->global_to_B,pcis->D,counter,ADD_VALUES,SCATTER_REVERSE);
256: VecScatterBegin(pcis->global_to_B,counter,pcis->vec1_B,INSERT_VALUES,SCATTER_FORWARD);
257: VecScatterEnd(pcis->global_to_B,counter,pcis->vec1_B,INSERT_VALUES,SCATTER_FORWARD);
258: VecPointwiseDivide(pcis->D,pcis->D,pcis->vec1_B);
259: VecDestroy(&counter);
261: /* See historical note 01, at the bottom of this file. */
263: /* Creating the KSP contexts for the local Dirichlet and Neumann problems */
264: if (computesolvers) {
265: PC pc_ctx;
267: pcis->pure_neumann = matis->pure_neumann;
268: /* Dirichlet */
269: KSPCreate(PETSC_COMM_SELF,&pcis->ksp_D);
270: KSPSetErrorIfNotConverged(pcis->ksp_D,pc->erroriffailure);
271: PetscObjectIncrementTabLevel((PetscObject)pcis->ksp_D,(PetscObject)pc,1);
272: KSPSetOperators(pcis->ksp_D,pcis->A_II,pcis->A_II);
273: KSPSetOptionsPrefix(pcis->ksp_D,"is_localD_");
274: KSPGetPC(pcis->ksp_D,&pc_ctx);
275: PCSetType(pc_ctx,PCLU);
276: KSPSetType(pcis->ksp_D,KSPPREONLY);
277: KSPSetFromOptions(pcis->ksp_D);
278: /* the vectors in the following line are dummy arguments, just telling the KSP the vector size. Values are not used */
279: KSPSetUp(pcis->ksp_D);
280: /* Neumann */
281: KSPCreate(PETSC_COMM_SELF,&pcis->ksp_N);
282: KSPSetErrorIfNotConverged(pcis->ksp_N,pc->erroriffailure);
283: PetscObjectIncrementTabLevel((PetscObject)pcis->ksp_N,(PetscObject)pc,1);
284: KSPSetOperators(pcis->ksp_N,matis->A,matis->A);
285: KSPSetOptionsPrefix(pcis->ksp_N,"is_localN_");
286: KSPGetPC(pcis->ksp_N,&pc_ctx);
287: PCSetType(pc_ctx,PCLU);
288: KSPSetType(pcis->ksp_N,KSPPREONLY);
289: KSPSetFromOptions(pcis->ksp_N);
290: {
291: PetscBool damp_fixed = PETSC_FALSE,
292: remove_nullspace_fixed = PETSC_FALSE,
293: set_damping_factor_floating = PETSC_FALSE,
294: not_damp_floating = PETSC_FALSE,
295: not_remove_nullspace_floating = PETSC_FALSE;
296: PetscReal fixed_factor,
297: floating_factor;
299: PetscOptionsGetReal(((PetscObject)pc_ctx)->options,((PetscObject)pc_ctx)->prefix,"-pc_is_damp_fixed",&fixed_factor,&damp_fixed);
300: if (!damp_fixed) fixed_factor = 0.0;
301: PetscOptionsGetBool(((PetscObject)pc_ctx)->options,((PetscObject)pc_ctx)->prefix,"-pc_is_damp_fixed",&damp_fixed,NULL);
303: PetscOptionsGetBool(((PetscObject)pc_ctx)->options,((PetscObject)pc_ctx)->prefix,"-pc_is_remove_nullspace_fixed",&remove_nullspace_fixed,NULL);
305: PetscOptionsGetReal(((PetscObject)pc_ctx)->options,((PetscObject)pc_ctx)->prefix,"-pc_is_set_damping_factor_floating",
306: &floating_factor,&set_damping_factor_floating);
307: if (!set_damping_factor_floating) floating_factor = 0.0;
308: PetscOptionsGetBool(((PetscObject)pc_ctx)->options,((PetscObject)pc_ctx)->prefix,"-pc_is_set_damping_factor_floating",&set_damping_factor_floating,NULL);
309: if (!set_damping_factor_floating) floating_factor = 1.e-12;
311: PetscOptionsGetBool(((PetscObject)pc_ctx)->options,((PetscObject)pc_ctx)->prefix,"-pc_is_not_damp_floating",¬_damp_floating,NULL);
313: PetscOptionsGetBool(((PetscObject)pc_ctx)->options,((PetscObject)pc_ctx)->prefix,"-pc_is_not_remove_nullspace_floating",¬_remove_nullspace_floating,NULL);
315: if (pcis->pure_neumann) { /* floating subdomain */
316: if (!(not_damp_floating)) {
317: PCFactorSetShiftType(pc_ctx,MAT_SHIFT_NONZERO);
318: PCFactorSetShiftAmount(pc_ctx,floating_factor);
319: }
320: if (!(not_remove_nullspace_floating)) {
321: MatNullSpace nullsp;
322: MatNullSpaceCreate(PETSC_COMM_SELF,PETSC_TRUE,0,NULL,&nullsp);
323: MatSetNullSpace(matis->A,nullsp);
324: MatNullSpaceDestroy(&nullsp);
325: }
326: } else { /* fixed subdomain */
327: if (damp_fixed) {
328: PCFactorSetShiftType(pc_ctx,MAT_SHIFT_NONZERO);
329: PCFactorSetShiftAmount(pc_ctx,floating_factor);
330: }
331: if (remove_nullspace_fixed) {
332: MatNullSpace nullsp;
333: MatNullSpaceCreate(PETSC_COMM_SELF,PETSC_TRUE,0,NULL,&nullsp);
334: MatSetNullSpace(matis->A,nullsp);
335: MatNullSpaceDestroy(&nullsp);
336: }
337: }
338: }
339: /* the vectors in the following line are dummy arguments, just telling the KSP the vector size. Values are not used */
340: KSPSetUp(pcis->ksp_N);
341: }
343: return(0);
344: }
346: /* -------------------------------------------------------------------------- */
347: /*
348: PCISDestroy -
349: */
350: PetscErrorCode PCISDestroy(PC pc)
351: {
352: PC_IS *pcis = (PC_IS*)(pc->data);
356: ISDestroy(&pcis->is_B_local);
357: ISDestroy(&pcis->is_I_local);
358: ISDestroy(&pcis->is_B_global);
359: ISDestroy(&pcis->is_I_global);
360: MatDestroy(&pcis->A_II);
361: MatDestroy(&pcis->A_IB);
362: MatDestroy(&pcis->A_BI);
363: MatDestroy(&pcis->A_BB);
364: VecDestroy(&pcis->D);
365: KSPDestroy(&pcis->ksp_N);
366: KSPDestroy(&pcis->ksp_D);
367: VecDestroy(&pcis->vec1_N);
368: VecDestroy(&pcis->vec2_N);
369: VecDestroy(&pcis->vec1_D);
370: VecDestroy(&pcis->vec2_D);
371: VecDestroy(&pcis->vec3_D);
372: VecDestroy(&pcis->vec4_D);
373: VecDestroy(&pcis->vec1_B);
374: VecDestroy(&pcis->vec2_B);
375: VecDestroy(&pcis->vec3_B);
376: VecDestroy(&pcis->vec1_global);
377: VecScatterDestroy(&pcis->global_to_D);
378: VecScatterDestroy(&pcis->N_to_B);
379: VecScatterDestroy(&pcis->N_to_D);
380: VecScatterDestroy(&pcis->global_to_B);
381: PetscFree(pcis->work_N);
382: if (pcis->n_neigh > -1) {
383: ISLocalToGlobalMappingRestoreInfo(pcis->mapping,&(pcis->n_neigh),&(pcis->neigh),&(pcis->n_shared),&(pcis->shared));
384: }
385: ISLocalToGlobalMappingDestroy(&pcis->mapping);
386: ISLocalToGlobalMappingDestroy(&pcis->BtoNmap);
387: PetscObjectComposeFunction((PetscObject)pc,"PCISSetUseStiffnessScaling_C",NULL);
388: PetscObjectComposeFunction((PetscObject)pc,"PCISSetSubdomainScalingFactor_C",NULL);
389: PetscObjectComposeFunction((PetscObject)pc,"PCISSetSubdomainDiagonalScaling_C",NULL);
390: return(0);
391: }
393: /* -------------------------------------------------------------------------- */
394: /*
395: PCISCreate -
396: */
397: PetscErrorCode PCISCreate(PC pc)
398: {
399: PC_IS *pcis = (PC_IS*)(pc->data);
403: pcis->is_B_local = 0;
404: pcis->is_I_local = 0;
405: pcis->is_B_global = 0;
406: pcis->is_I_global = 0;
407: pcis->A_II = 0;
408: pcis->A_IB = 0;
409: pcis->A_BI = 0;
410: pcis->A_BB = 0;
411: pcis->D = 0;
412: pcis->ksp_N = 0;
413: pcis->ksp_D = 0;
414: pcis->vec1_N = 0;
415: pcis->vec2_N = 0;
416: pcis->vec1_D = 0;
417: pcis->vec2_D = 0;
418: pcis->vec3_D = 0;
419: pcis->vec1_B = 0;
420: pcis->vec2_B = 0;
421: pcis->vec3_B = 0;
422: pcis->vec1_global = 0;
423: pcis->work_N = 0;
424: pcis->global_to_D = 0;
425: pcis->N_to_B = 0;
426: pcis->N_to_D = 0;
427: pcis->global_to_B = 0;
428: pcis->mapping = 0;
429: pcis->BtoNmap = 0;
430: pcis->n_neigh = -1;
431: pcis->scaling_factor = 1.0;
432: pcis->reusesubmatrices = PETSC_TRUE;
433: /* composing functions */
434: PetscObjectComposeFunction((PetscObject)pc,"PCISSetUseStiffnessScaling_C",PCISSetUseStiffnessScaling_IS);
435: PetscObjectComposeFunction((PetscObject)pc,"PCISSetSubdomainScalingFactor_C",PCISSetSubdomainScalingFactor_IS);
436: PetscObjectComposeFunction((PetscObject)pc,"PCISSetSubdomainDiagonalScaling_C",PCISSetSubdomainDiagonalScaling_IS);
437: return(0);
438: }
440: /* -------------------------------------------------------------------------- */
441: /*
442: PCISApplySchur -
444: Input parameters:
445: . pc - preconditioner context
446: . v - vector to which the Schur complement is to be applied (it is NOT modified inside this function, UNLESS vec2_B is null)
448: Output parameters:
449: . vec1_B - result of Schur complement applied to chunk
450: . vec2_B - garbage (used as work space), or null (and v is used as workspace)
451: . vec1_D - garbage (used as work space)
452: . vec2_D - garbage (used as work space)
454: */
455: PetscErrorCode PCISApplySchur(PC pc, Vec v, Vec vec1_B, Vec vec2_B, Vec vec1_D, Vec vec2_D)
456: {
458: PC_IS *pcis = (PC_IS*)(pc->data);
461: if (!vec2_B) vec2_B = v;
463: MatMult(pcis->A_BB,v,vec1_B);
464: MatMult(pcis->A_IB,v,vec1_D);
465: KSPSolve(pcis->ksp_D,vec1_D,vec2_D);
466: MatMult(pcis->A_BI,vec2_D,vec2_B);
467: VecAXPY(vec1_B,-1.0,vec2_B);
468: return(0);
469: }
471: /* -------------------------------------------------------------------------- */
472: /*
473: PCISScatterArrayNToVecB - Scatters interface node values from a big array (of all local nodes, interior or interface,
474: including ghosts) into an interface vector, when in SCATTER_FORWARD mode, or vice-versa, when in SCATTER_REVERSE
475: mode.
477: Input parameters:
478: . pc - preconditioner context
479: . array_N - [when in SCATTER_FORWARD mode] Array to be scattered into the vector
480: . v_B - [when in SCATTER_REVERSE mode] Vector to be scattered into the array
482: Output parameter:
483: . array_N - [when in SCATTER_REVERSE mode] Array to receive the scattered vector
484: . v_B - [when in SCATTER_FORWARD mode] Vector to receive the scattered array
486: Notes:
487: The entries in the array that do not correspond to interface nodes remain unaltered.
488: */
489: PetscErrorCode PCISScatterArrayNToVecB(PetscScalar *array_N, Vec v_B, InsertMode imode, ScatterMode smode, PC pc)
490: {
491: PetscInt i;
492: const PetscInt *idex;
494: PetscScalar *array_B;
495: PC_IS *pcis = (PC_IS*)(pc->data);
498: VecGetArray(v_B,&array_B);
499: ISGetIndices(pcis->is_B_local,&idex);
501: if (smode == SCATTER_FORWARD) {
502: if (imode == INSERT_VALUES) {
503: for (i=0; i<pcis->n_B; i++) array_B[i] = array_N[idex[i]];
504: } else { /* ADD_VALUES */
505: for (i=0; i<pcis->n_B; i++) array_B[i] += array_N[idex[i]];
506: }
507: } else { /* SCATTER_REVERSE */
508: if (imode == INSERT_VALUES) {
509: for (i=0; i<pcis->n_B; i++) array_N[idex[i]] = array_B[i];
510: } else { /* ADD_VALUES */
511: for (i=0; i<pcis->n_B; i++) array_N[idex[i]] += array_B[i];
512: }
513: }
514: ISRestoreIndices(pcis->is_B_local,&idex);
515: VecRestoreArray(v_B,&array_B);
516: return(0);
517: }
519: /* -------------------------------------------------------------------------- */
520: /*
521: PCISApplyInvSchur - Solves the Neumann problem related to applying the inverse of the Schur complement.
522: More precisely, solves the problem:
523: [ A_II A_IB ] [ . ] [ 0 ]
524: [ ] [ ] = [ ]
525: [ A_BI A_BB ] [ x ] [ b ]
527: Input parameters:
528: . pc - preconditioner context
529: . b - vector of local interface nodes (including ghosts)
531: Output parameters:
532: . x - vector of local interface nodes (including ghosts); returns the application of the inverse of the Schur
533: complement to b
534: . vec1_N - vector of local nodes (interior and interface, including ghosts); returns garbage (used as work space)
535: . vec2_N - vector of local nodes (interior and interface, including ghosts); returns garbage (used as work space)
537: */
538: PetscErrorCode PCISApplyInvSchur(PC pc, Vec b, Vec x, Vec vec1_N, Vec vec2_N)
539: {
541: PC_IS *pcis = (PC_IS*)(pc->data);
544: /*
545: Neumann solvers.
546: Applying the inverse of the local Schur complement, i.e, solving a Neumann
547: Problem with zero at the interior nodes of the RHS and extracting the interface
548: part of the solution. inverse Schur complement is applied to b and the result
549: is stored in x.
550: */
551: /* Setting the RHS vec1_N */
552: VecSet(vec1_N,0.0);
553: VecScatterBegin(pcis->N_to_B,b,vec1_N,INSERT_VALUES,SCATTER_REVERSE);
554: VecScatterEnd (pcis->N_to_B,b,vec1_N,INSERT_VALUES,SCATTER_REVERSE);
555: /* Checking for consistency of the RHS */
556: {
557: PetscBool flg = PETSC_FALSE;
558: PetscOptionsGetBool(NULL,NULL,"-pc_is_check_consistency",&flg,NULL);
559: if (flg) {
560: PetscScalar average;
561: PetscViewer viewer;
562: PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)pc),&viewer);
564: VecSum(vec1_N,&average);
565: average = average / ((PetscReal)pcis->n);
566: PetscViewerASCIIPushSynchronized(viewer);
567: if (pcis->pure_neumann) {
568: PetscViewerASCIISynchronizedPrintf(viewer,"Subdomain %04d is floating. Average = % 1.14e\n",PetscGlobalRank,PetscAbsScalar(average));
569: } else {
570: PetscViewerASCIISynchronizedPrintf(viewer,"Subdomain %04d is fixed. Average = % 1.14e\n",PetscGlobalRank,PetscAbsScalar(average));
571: }
572: PetscViewerFlush(viewer);
573: PetscViewerASCIIPopSynchronized(viewer);
574: }
575: }
576: /* Solving the system for vec2_N */
577: KSPSolve(pcis->ksp_N,vec1_N,vec2_N);
578: /* Extracting the local interface vector out of the solution */
579: VecScatterBegin(pcis->N_to_B,vec2_N,x,INSERT_VALUES,SCATTER_FORWARD);
580: VecScatterEnd (pcis->N_to_B,vec2_N,x,INSERT_VALUES,SCATTER_FORWARD);
581: return(0);
582: }