Actual source code: pcis.c
petsc-3.10.5 2019-03-28
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 computematrices, PetscBool computesolvers)
118: {
119: PC_IS *pcis = (PC_IS*)(pc->data);
120: Mat_IS *matis;
121: MatReuse reuse;
123: PetscBool flg,issbaij;
126: PetscObjectTypeCompare((PetscObject)pc->pmat,MATIS,&flg);
127: if (!flg) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONG,"Preconditioner type of Neumann Neumman requires matrix of type MATIS");
128: matis = (Mat_IS*)pc->pmat->data;
130: /* first time creation, get info on substructuring */
131: if (!pc->setupcalled) {
132: PetscInt n_I;
133: PetscInt *idx_I_local,*idx_B_local,*idx_I_global,*idx_B_global;
134: PetscBT bt;
135: PetscInt i,j;
137: /* get info on mapping */
138: PetscObjectReference((PetscObject)pc->pmat->rmap->mapping);
139: ISLocalToGlobalMappingDestroy(&pcis->mapping);
140: pcis->mapping = pc->pmat->rmap->mapping;
141: ISLocalToGlobalMappingGetSize(pcis->mapping,&pcis->n);
142: ISLocalToGlobalMappingGetInfo(pcis->mapping,&(pcis->n_neigh),&(pcis->neigh),&(pcis->n_shared),&(pcis->shared));
144: /* Identifying interior and interface nodes, in local numbering */
145: PetscBTCreate(pcis->n,&bt);
146: for (i=0;i<pcis->n_neigh;i++)
147: for (j=0;j<pcis->n_shared[i];j++) {
148: PetscBTSet(bt,pcis->shared[i][j]);
149: }
151: /* Creating local and global index sets for interior and inteface nodes. */
152: PetscMalloc1(pcis->n,&idx_I_local);
153: PetscMalloc1(pcis->n,&idx_B_local);
154: for (i=0, pcis->n_B=0, n_I=0; i<pcis->n; i++) {
155: if (!PetscBTLookup(bt,i)) {
156: idx_I_local[n_I] = i;
157: n_I++;
158: } else {
159: idx_B_local[pcis->n_B] = i;
160: pcis->n_B++;
161: }
162: }
164: /* Getting the global numbering */
165: idx_B_global = idx_I_local + n_I; /* Just avoiding allocating extra memory, since we have vacant space */
166: idx_I_global = idx_B_local + pcis->n_B;
167: ISLocalToGlobalMappingApply(pcis->mapping,pcis->n_B,idx_B_local,idx_B_global);
168: ISLocalToGlobalMappingApply(pcis->mapping,n_I,idx_I_local,idx_I_global);
170: /* Creating the index sets */
171: ISCreateGeneral(PETSC_COMM_SELF,pcis->n_B,idx_B_local,PETSC_COPY_VALUES, &pcis->is_B_local);
172: ISCreateGeneral(PetscObjectComm((PetscObject)pc),pcis->n_B,idx_B_global,PETSC_COPY_VALUES,&pcis->is_B_global);
173: ISCreateGeneral(PETSC_COMM_SELF,n_I,idx_I_local,PETSC_COPY_VALUES, &pcis->is_I_local);
174: ISCreateGeneral(PetscObjectComm((PetscObject)pc),n_I,idx_I_global,PETSC_COPY_VALUES,&pcis->is_I_global);
176: /* Freeing memory */
177: PetscFree(idx_B_local);
178: PetscFree(idx_I_local);
179: PetscBTDestroy(&bt);
181: /* Creating work vectors and arrays */
182: VecDuplicate(matis->x,&pcis->vec1_N);
183: VecDuplicate(pcis->vec1_N,&pcis->vec2_N);
184: VecCreate(PETSC_COMM_SELF,&pcis->vec1_D);
185: VecSetSizes(pcis->vec1_D,pcis->n-pcis->n_B,PETSC_DECIDE);
186: VecSetType(pcis->vec1_D,((PetscObject)pcis->vec1_N)->type_name);
187: VecDuplicate(pcis->vec1_D,&pcis->vec2_D);
188: VecDuplicate(pcis->vec1_D,&pcis->vec3_D);
189: VecDuplicate(pcis->vec1_D,&pcis->vec4_D);
190: VecCreate(PETSC_COMM_SELF,&pcis->vec1_B);
191: VecSetSizes(pcis->vec1_B,pcis->n_B,PETSC_DECIDE);
192: VecSetType(pcis->vec1_B,((PetscObject)pcis->vec1_N)->type_name);
193: VecDuplicate(pcis->vec1_B,&pcis->vec2_B);
194: VecDuplicate(pcis->vec1_B,&pcis->vec3_B);
195: MatCreateVecs(pc->pmat,&pcis->vec1_global,0);
196: PetscMalloc1(pcis->n,&pcis->work_N);
197: /* scaling vector */
198: if (!pcis->D) { /* it can happen that the user passed in a scaling vector via PCISSetSubdomainDiagonalScaling */
199: VecDuplicate(pcis->vec1_B,&pcis->D);
200: VecSet(pcis->D,pcis->scaling_factor);
201: }
203: /* Creating the scatter contexts */
204: VecScatterCreate(pcis->vec1_N,pcis->is_I_local,pcis->vec1_D,(IS)0,&pcis->N_to_D);
205: VecScatterCreate(pcis->vec1_global,pcis->is_I_global,pcis->vec1_D,(IS)0,&pcis->global_to_D);
206: VecScatterCreate(pcis->vec1_N,pcis->is_B_local,pcis->vec1_B,(IS)0,&pcis->N_to_B);
207: VecScatterCreate(pcis->vec1_global,pcis->is_B_global,pcis->vec1_B,(IS)0,&pcis->global_to_B);
209: /* map from boundary to local */
210: ISLocalToGlobalMappingCreateIS(pcis->is_B_local,&pcis->BtoNmap);
211: }
213: /*
214: Extracting the blocks A_II, A_BI, A_IB and A_BB from A. If the numbering
215: is such that interior nodes come first than the interface ones, we have
217: [ A_II | A_IB ]
218: A = [------+------]
219: [ A_BI | A_BB ]
220: */
221: if (computematrices) {
222: reuse = MAT_INITIAL_MATRIX;
223: if (pcis->reusesubmatrices && pc->setupcalled) {
224: if (pc->flag == SAME_NONZERO_PATTERN) {
225: reuse = MAT_REUSE_MATRIX;
226: } else {
227: reuse = MAT_INITIAL_MATRIX;
228: }
229: }
230: if (reuse == MAT_INITIAL_MATRIX) {
231: MatDestroy(&pcis->A_II);
232: MatDestroy(&pcis->A_IB);
233: MatDestroy(&pcis->A_BI);
234: MatDestroy(&pcis->A_BB);
235: }
237: MatCreateSubMatrix(matis->A,pcis->is_I_local,pcis->is_I_local,reuse,&pcis->A_II);
238: MatCreateSubMatrix(matis->A,pcis->is_B_local,pcis->is_B_local,reuse,&pcis->A_BB);
239: PetscObjectTypeCompare((PetscObject)matis->A,MATSEQSBAIJ,&issbaij);
240: if (!issbaij) {
241: MatCreateSubMatrix(matis->A,pcis->is_I_local,pcis->is_B_local,reuse,&pcis->A_IB);
242: MatCreateSubMatrix(matis->A,pcis->is_B_local,pcis->is_I_local,reuse,&pcis->A_BI);
243: } else {
244: Mat newmat;
246: MatConvert(matis->A,MATSEQBAIJ,MAT_INITIAL_MATRIX,&newmat);
247: MatCreateSubMatrix(newmat,pcis->is_I_local,pcis->is_B_local,reuse,&pcis->A_IB);
248: MatCreateSubMatrix(newmat,pcis->is_B_local,pcis->is_I_local,reuse,&pcis->A_BI);
249: MatDestroy(&newmat);
250: }
251: }
253: /* Creating scaling vector D */
254: PetscOptionsGetBool(((PetscObject)pc)->options,((PetscObject)pc)->prefix,"-pc_is_use_stiffness_scaling",&pcis->use_stiffness_scaling,NULL);
255: if (pcis->use_stiffness_scaling) {
256: PetscScalar *a;
257: PetscInt i,n;
259: if (pcis->A_BB) {
260: MatGetDiagonal(pcis->A_BB,pcis->D);
261: } else {
262: MatGetDiagonal(matis->A,pcis->vec1_N);
263: VecScatterBegin(pcis->N_to_B,pcis->vec1_N,pcis->D,INSERT_VALUES,SCATTER_FORWARD);
264: VecScatterEnd(pcis->N_to_B,pcis->vec1_N,pcis->D,INSERT_VALUES,SCATTER_FORWARD);
265: }
266: VecAbs(pcis->D);
267: VecGetLocalSize(pcis->D,&n);
268: VecGetArray(pcis->D,&a);
269: for (i=0;i<n;i++) if (PetscAbsScalar(a[i])<PETSC_SMALL) a[i] = 1.0;
270: VecRestoreArray(pcis->D,&a);
271: }
272: VecSet(pcis->vec1_global,0.0);
273: VecScatterBegin(pcis->global_to_B,pcis->D,pcis->vec1_global,ADD_VALUES,SCATTER_REVERSE);
274: VecScatterEnd(pcis->global_to_B,pcis->D,pcis->vec1_global,ADD_VALUES,SCATTER_REVERSE);
275: VecScatterBegin(pcis->global_to_B,pcis->vec1_global,pcis->vec1_B,INSERT_VALUES,SCATTER_FORWARD);
276: VecScatterEnd(pcis->global_to_B,pcis->vec1_global,pcis->vec1_B,INSERT_VALUES,SCATTER_FORWARD);
277: VecPointwiseDivide(pcis->D,pcis->D,pcis->vec1_B);
279: /* See historical note 01, at the bottom of this file. */
281: /* Creating the KSP contexts for the local Dirichlet and Neumann problems */
282: if (computesolvers) {
283: PC pc_ctx;
285: pcis->pure_neumann = matis->pure_neumann;
286: /* Dirichlet */
287: KSPCreate(PETSC_COMM_SELF,&pcis->ksp_D);
288: KSPSetErrorIfNotConverged(pcis->ksp_D,pc->erroriffailure);
289: PetscObjectIncrementTabLevel((PetscObject)pcis->ksp_D,(PetscObject)pc,1);
290: KSPSetOperators(pcis->ksp_D,pcis->A_II,pcis->A_II);
291: KSPSetOptionsPrefix(pcis->ksp_D,"is_localD_");
292: KSPGetPC(pcis->ksp_D,&pc_ctx);
293: PCSetType(pc_ctx,PCLU);
294: KSPSetType(pcis->ksp_D,KSPPREONLY);
295: KSPSetFromOptions(pcis->ksp_D);
296: /* the vectors in the following line are dummy arguments, just telling the KSP the vector size. Values are not used */
297: KSPSetUp(pcis->ksp_D);
298: /* Neumann */
299: KSPCreate(PETSC_COMM_SELF,&pcis->ksp_N);
300: KSPSetErrorIfNotConverged(pcis->ksp_N,pc->erroriffailure);
301: PetscObjectIncrementTabLevel((PetscObject)pcis->ksp_N,(PetscObject)pc,1);
302: KSPSetOperators(pcis->ksp_N,matis->A,matis->A);
303: KSPSetOptionsPrefix(pcis->ksp_N,"is_localN_");
304: KSPGetPC(pcis->ksp_N,&pc_ctx);
305: PCSetType(pc_ctx,PCLU);
306: KSPSetType(pcis->ksp_N,KSPPREONLY);
307: KSPSetFromOptions(pcis->ksp_N);
308: {
309: PetscBool damp_fixed = PETSC_FALSE,
310: remove_nullspace_fixed = PETSC_FALSE,
311: set_damping_factor_floating = PETSC_FALSE,
312: not_damp_floating = PETSC_FALSE,
313: not_remove_nullspace_floating = PETSC_FALSE;
314: PetscReal fixed_factor,
315: floating_factor;
317: PetscOptionsGetReal(((PetscObject)pc_ctx)->options,((PetscObject)pc_ctx)->prefix,"-pc_is_damp_fixed",&fixed_factor,&damp_fixed);
318: if (!damp_fixed) fixed_factor = 0.0;
319: PetscOptionsGetBool(((PetscObject)pc_ctx)->options,((PetscObject)pc_ctx)->prefix,"-pc_is_damp_fixed",&damp_fixed,NULL);
321: PetscOptionsGetBool(((PetscObject)pc_ctx)->options,((PetscObject)pc_ctx)->prefix,"-pc_is_remove_nullspace_fixed",&remove_nullspace_fixed,NULL);
323: PetscOptionsGetReal(((PetscObject)pc_ctx)->options,((PetscObject)pc_ctx)->prefix,"-pc_is_set_damping_factor_floating",
324: &floating_factor,&set_damping_factor_floating);
325: if (!set_damping_factor_floating) floating_factor = 0.0;
326: PetscOptionsGetBool(((PetscObject)pc_ctx)->options,((PetscObject)pc_ctx)->prefix,"-pc_is_set_damping_factor_floating",&set_damping_factor_floating,NULL);
327: if (!set_damping_factor_floating) floating_factor = 1.e-12;
329: PetscOptionsGetBool(((PetscObject)pc_ctx)->options,((PetscObject)pc_ctx)->prefix,"-pc_is_not_damp_floating",¬_damp_floating,NULL);
331: PetscOptionsGetBool(((PetscObject)pc_ctx)->options,((PetscObject)pc_ctx)->prefix,"-pc_is_not_remove_nullspace_floating",¬_remove_nullspace_floating,NULL);
333: if (pcis->pure_neumann) { /* floating subdomain */
334: if (!(not_damp_floating)) {
335: PCFactorSetShiftType(pc_ctx,MAT_SHIFT_NONZERO);
336: PCFactorSetShiftAmount(pc_ctx,floating_factor);
337: }
338: if (!(not_remove_nullspace_floating)) {
339: MatNullSpace nullsp;
340: MatNullSpaceCreate(PETSC_COMM_SELF,PETSC_TRUE,0,NULL,&nullsp);
341: MatSetNullSpace(matis->A,nullsp);
342: MatNullSpaceDestroy(&nullsp);
343: }
344: } else { /* fixed subdomain */
345: if (damp_fixed) {
346: PCFactorSetShiftType(pc_ctx,MAT_SHIFT_NONZERO);
347: PCFactorSetShiftAmount(pc_ctx,floating_factor);
348: }
349: if (remove_nullspace_fixed) {
350: MatNullSpace nullsp;
351: MatNullSpaceCreate(PETSC_COMM_SELF,PETSC_TRUE,0,NULL,&nullsp);
352: MatSetNullSpace(matis->A,nullsp);
353: MatNullSpaceDestroy(&nullsp);
354: }
355: }
356: }
357: /* the vectors in the following line are dummy arguments, just telling the KSP the vector size. Values are not used */
358: KSPSetUp(pcis->ksp_N);
359: }
360: return(0);
361: }
363: /* -------------------------------------------------------------------------- */
364: /*
365: PCISDestroy -
366: */
367: PetscErrorCode PCISDestroy(PC pc)
368: {
369: PC_IS *pcis = (PC_IS*)(pc->data);
373: ISDestroy(&pcis->is_B_local);
374: ISDestroy(&pcis->is_I_local);
375: ISDestroy(&pcis->is_B_global);
376: ISDestroy(&pcis->is_I_global);
377: MatDestroy(&pcis->A_II);
378: MatDestroy(&pcis->A_IB);
379: MatDestroy(&pcis->A_BI);
380: MatDestroy(&pcis->A_BB);
381: VecDestroy(&pcis->D);
382: KSPDestroy(&pcis->ksp_N);
383: KSPDestroy(&pcis->ksp_D);
384: VecDestroy(&pcis->vec1_N);
385: VecDestroy(&pcis->vec2_N);
386: VecDestroy(&pcis->vec1_D);
387: VecDestroy(&pcis->vec2_D);
388: VecDestroy(&pcis->vec3_D);
389: VecDestroy(&pcis->vec4_D);
390: VecDestroy(&pcis->vec1_B);
391: VecDestroy(&pcis->vec2_B);
392: VecDestroy(&pcis->vec3_B);
393: VecDestroy(&pcis->vec1_global);
394: VecScatterDestroy(&pcis->global_to_D);
395: VecScatterDestroy(&pcis->N_to_B);
396: VecScatterDestroy(&pcis->N_to_D);
397: VecScatterDestroy(&pcis->global_to_B);
398: PetscFree(pcis->work_N);
399: if (pcis->n_neigh > -1) {
400: ISLocalToGlobalMappingRestoreInfo(pcis->mapping,&(pcis->n_neigh),&(pcis->neigh),&(pcis->n_shared),&(pcis->shared));
401: }
402: ISLocalToGlobalMappingDestroy(&pcis->mapping);
403: ISLocalToGlobalMappingDestroy(&pcis->BtoNmap);
404: PetscObjectComposeFunction((PetscObject)pc,"PCISSetUseStiffnessScaling_C",NULL);
405: PetscObjectComposeFunction((PetscObject)pc,"PCISSetSubdomainScalingFactor_C",NULL);
406: PetscObjectComposeFunction((PetscObject)pc,"PCISSetSubdomainDiagonalScaling_C",NULL);
407: return(0);
408: }
410: /* -------------------------------------------------------------------------- */
411: /*
412: PCISCreate -
413: */
414: PetscErrorCode PCISCreate(PC pc)
415: {
416: PC_IS *pcis = (PC_IS*)(pc->data);
420: pcis->n_neigh = -1;
421: pcis->scaling_factor = 1.0;
422: pcis->reusesubmatrices = PETSC_TRUE;
423: /* composing functions */
424: PetscObjectComposeFunction((PetscObject)pc,"PCISSetUseStiffnessScaling_C",PCISSetUseStiffnessScaling_IS);
425: PetscObjectComposeFunction((PetscObject)pc,"PCISSetSubdomainScalingFactor_C",PCISSetSubdomainScalingFactor_IS);
426: PetscObjectComposeFunction((PetscObject)pc,"PCISSetSubdomainDiagonalScaling_C",PCISSetSubdomainDiagonalScaling_IS);
427: return(0);
428: }
430: /* -------------------------------------------------------------------------- */
431: /*
432: PCISApplySchur -
434: Input parameters:
435: . pc - preconditioner context
436: . v - vector to which the Schur complement is to be applied (it is NOT modified inside this function, UNLESS vec2_B is null)
438: Output parameters:
439: . vec1_B - result of Schur complement applied to chunk
440: . vec2_B - garbage (used as work space), or null (and v is used as workspace)
441: . vec1_D - garbage (used as work space)
442: . vec2_D - garbage (used as work space)
444: */
445: PetscErrorCode PCISApplySchur(PC pc, Vec v, Vec vec1_B, Vec vec2_B, Vec vec1_D, Vec vec2_D)
446: {
448: PC_IS *pcis = (PC_IS*)(pc->data);
451: if (!vec2_B) vec2_B = v;
453: MatMult(pcis->A_BB,v,vec1_B);
454: MatMult(pcis->A_IB,v,vec1_D);
455: KSPSolve(pcis->ksp_D,vec1_D,vec2_D);
456: MatMult(pcis->A_BI,vec2_D,vec2_B);
457: VecAXPY(vec1_B,-1.0,vec2_B);
458: return(0);
459: }
461: /* -------------------------------------------------------------------------- */
462: /*
463: PCISScatterArrayNToVecB - Scatters interface node values from a big array (of all local nodes, interior or interface,
464: including ghosts) into an interface vector, when in SCATTER_FORWARD mode, or vice-versa, when in SCATTER_REVERSE
465: mode.
467: Input parameters:
468: . pc - preconditioner context
469: . array_N - [when in SCATTER_FORWARD mode] Array to be scattered into the vector
470: . v_B - [when in SCATTER_REVERSE mode] Vector to be scattered into the array
472: Output parameter:
473: . array_N - [when in SCATTER_REVERSE mode] Array to receive the scattered vector
474: . v_B - [when in SCATTER_FORWARD mode] Vector to receive the scattered array
476: Notes:
477: The entries in the array that do not correspond to interface nodes remain unaltered.
478: */
479: PetscErrorCode PCISScatterArrayNToVecB(PetscScalar *array_N, Vec v_B, InsertMode imode, ScatterMode smode, PC pc)
480: {
481: PetscInt i;
482: const PetscInt *idex;
484: PetscScalar *array_B;
485: PC_IS *pcis = (PC_IS*)(pc->data);
488: VecGetArray(v_B,&array_B);
489: ISGetIndices(pcis->is_B_local,&idex);
491: if (smode == SCATTER_FORWARD) {
492: if (imode == INSERT_VALUES) {
493: for (i=0; i<pcis->n_B; i++) array_B[i] = array_N[idex[i]];
494: } else { /* ADD_VALUES */
495: for (i=0; i<pcis->n_B; i++) array_B[i] += array_N[idex[i]];
496: }
497: } else { /* SCATTER_REVERSE */
498: if (imode == INSERT_VALUES) {
499: for (i=0; i<pcis->n_B; i++) array_N[idex[i]] = array_B[i];
500: } else { /* ADD_VALUES */
501: for (i=0; i<pcis->n_B; i++) array_N[idex[i]] += array_B[i];
502: }
503: }
504: ISRestoreIndices(pcis->is_B_local,&idex);
505: VecRestoreArray(v_B,&array_B);
506: return(0);
507: }
509: /* -------------------------------------------------------------------------- */
510: /*
511: PCISApplyInvSchur - Solves the Neumann problem related to applying the inverse of the Schur complement.
512: More precisely, solves the problem:
513: [ A_II A_IB ] [ . ] [ 0 ]
514: [ ] [ ] = [ ]
515: [ A_BI A_BB ] [ x ] [ b ]
517: Input parameters:
518: . pc - preconditioner context
519: . b - vector of local interface nodes (including ghosts)
521: Output parameters:
522: . x - vector of local interface nodes (including ghosts); returns the application of the inverse of the Schur
523: complement to b
524: . vec1_N - vector of local nodes (interior and interface, including ghosts); returns garbage (used as work space)
525: . vec2_N - vector of local nodes (interior and interface, including ghosts); returns garbage (used as work space)
527: */
528: PetscErrorCode PCISApplyInvSchur(PC pc, Vec b, Vec x, Vec vec1_N, Vec vec2_N)
529: {
531: PC_IS *pcis = (PC_IS*)(pc->data);
534: /*
535: Neumann solvers.
536: Applying the inverse of the local Schur complement, i.e, solving a Neumann
537: Problem with zero at the interior nodes of the RHS and extracting the interface
538: part of the solution. inverse Schur complement is applied to b and the result
539: is stored in x.
540: */
541: /* Setting the RHS vec1_N */
542: VecSet(vec1_N,0.0);
543: VecScatterBegin(pcis->N_to_B,b,vec1_N,INSERT_VALUES,SCATTER_REVERSE);
544: VecScatterEnd (pcis->N_to_B,b,vec1_N,INSERT_VALUES,SCATTER_REVERSE);
545: /* Checking for consistency of the RHS */
546: {
547: PetscBool flg = PETSC_FALSE;
548: PetscOptionsGetBool(NULL,NULL,"-pc_is_check_consistency",&flg,NULL);
549: if (flg) {
550: PetscScalar average;
551: PetscViewer viewer;
552: PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)pc),&viewer);
554: VecSum(vec1_N,&average);
555: average = average / ((PetscReal)pcis->n);
556: PetscViewerASCIIPushSynchronized(viewer);
557: if (pcis->pure_neumann) {
558: PetscViewerASCIISynchronizedPrintf(viewer,"Subdomain %04d is floating. Average = % 1.14e\n",PetscGlobalRank,PetscAbsScalar(average));
559: } else {
560: PetscViewerASCIISynchronizedPrintf(viewer,"Subdomain %04d is fixed. Average = % 1.14e\n",PetscGlobalRank,PetscAbsScalar(average));
561: }
562: PetscViewerFlush(viewer);
563: PetscViewerASCIIPopSynchronized(viewer);
564: }
565: }
566: /* Solving the system for vec2_N */
567: KSPSolve(pcis->ksp_N,vec1_N,vec2_N);
568: /* Extracting the local interface vector out of the solution */
569: VecScatterBegin(pcis->N_to_B,vec2_N,x,INSERT_VALUES,SCATTER_FORWARD);
570: VecScatterEnd (pcis->N_to_B,vec2_N,x,INSERT_VALUES,SCATTER_FORWARD);
571: return(0);
572: }