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", ¬_damp_floating, NULL);
359: PetscOptionsGetBool(((PetscObject)pc_ctx)->options, ((PetscObject)pc_ctx)->prefix, "-pc_is_not_remove_nullspace_floating", ¬_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: }