Actual source code: pcis.c
1: #include <petsc/private/pcisimpl.h>
3: static PetscErrorCode PCISSetUseStiffnessScaling_IS(PC pc, PetscBool use)
4: {
5: PC_IS *pcis = (PC_IS *)pc->data;
7: PetscFunctionBegin;
8: pcis->use_stiffness_scaling = use;
9: PetscFunctionReturn(PETSC_SUCCESS);
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 it should use matrix diagonal to build partition of unity.
22: Level: intermediate
24: .seealso: [](ch_ksp), `PCBDDC`, `PCNN`, `PCISSetSubdomainDiagonalScaling()`, `PCISScatterArrayNToVecB()`,
25: `PCISSetSubdomainScalingFactor()`,
26: `PCISReset()`, `PCISInitialize()`, `PCISApplyInvSchur()`, `PCISApplySchur()`
27: @*/
28: PetscErrorCode PCISSetUseStiffnessScaling(PC pc, PetscBool use)
29: {
30: PetscFunctionBegin;
33: PetscTryMethod(pc, "PCISSetUseStiffnessScaling_C", (PC, PetscBool), (pc, use));
34: PetscFunctionReturn(PETSC_SUCCESS);
35: }
37: static PetscErrorCode PCISSetSubdomainDiagonalScaling_IS(PC pc, Vec scaling_factors)
38: {
39: PC_IS *pcis = (PC_IS *)pc->data;
41: PetscFunctionBegin;
42: PetscCall(PetscObjectReference((PetscObject)scaling_factors));
43: PetscCall(VecDestroy(&pcis->D));
44: pcis->D = scaling_factors;
45: if (pc->setupcalled) {
46: PetscInt sn;
48: PetscCall(VecGetSize(pcis->D, &sn));
49: if (sn == pcis->n) {
50: PetscCall(VecScatterBegin(pcis->N_to_B, pcis->D, pcis->vec1_B, INSERT_VALUES, SCATTER_FORWARD));
51: PetscCall(VecScatterEnd(pcis->N_to_B, pcis->D, pcis->vec1_B, INSERT_VALUES, SCATTER_FORWARD));
52: PetscCall(VecDestroy(&pcis->D));
53: PetscCall(VecDuplicate(pcis->vec1_B, &pcis->D));
54: PetscCall(VecCopy(pcis->vec1_B, pcis->D));
55: } else PetscCheck(sn == pcis->n_B, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Invalid size for scaling vector. Expected %" PetscInt_FMT " (or full %" PetscInt_FMT "), found %" PetscInt_FMT, pcis->n_B, pcis->n, sn);
56: }
57: PetscFunctionReturn(PETSC_SUCCESS);
58: }
60: /*@
61: PCISSetSubdomainDiagonalScaling - Set diagonal scaling for `PCIS`.
63: Logically Collective
65: Input Parameters:
66: + pc - the preconditioning context
67: - scaling_factors - scaling factors for the subdomain
69: Level: intermediate
71: Note:
72: Intended for use with jumping coefficients cases.
74: .seealso: [](ch_ksp), `PCBDDC`, `PCNN`, `PCISScatterArrayNToVecB()`,
75: `PCISSetSubdomainScalingFactor()`, `PCISSetUseStiffnessScaling()`,
76: `PCISReset()`, `PCISInitialize()`, `PCISApplyInvSchur()`, `PCISApplySchur()`
77: @*/
78: PetscErrorCode PCISSetSubdomainDiagonalScaling(PC pc, Vec scaling_factors)
79: {
80: PetscFunctionBegin;
83: PetscTryMethod(pc, "PCISSetSubdomainDiagonalScaling_C", (PC, Vec), (pc, scaling_factors));
84: PetscFunctionReturn(PETSC_SUCCESS);
85: }
87: static PetscErrorCode PCISSetSubdomainScalingFactor_IS(PC pc, PetscScalar scal)
88: {
89: PC_IS *pcis = (PC_IS *)pc->data;
91: PetscFunctionBegin;
92: pcis->scaling_factor = scal;
93: if (pcis->D) PetscCall(VecSet(pcis->D, pcis->scaling_factor));
94: PetscFunctionReturn(PETSC_SUCCESS);
95: }
97: /*@
98: PCISSetSubdomainScalingFactor - Set scaling factor for `PCIS`.
100: Not Collective
102: Input Parameters:
103: + pc - the preconditioning context
104: - scal - scaling factor for the subdomain
106: Level: intermediate
108: Note:
109: Intended for use with the jumping coefficients cases.
111: .seealso: [](ch_ksp), `PCBDDC`, `PCNN`, `PCISScatterArrayNToVecB()`,
112: `PCISSetSubdomainDiagonalScaling()`, `PCISSetUseStiffnessScaling()`,
113: `PCISReset()`, `PCISInitialize()`, `PCISApplyInvSchur()`, `PCISApplySchur()`
114: @*/
115: PetscErrorCode PCISSetSubdomainScalingFactor(PC pc, PetscScalar scal)
116: {
117: PetscFunctionBegin;
119: PetscTryMethod(pc, "PCISSetSubdomainScalingFactor_C", (PC, PetscScalar), (pc, scal));
120: PetscFunctionReturn(PETSC_SUCCESS);
121: }
123: /*@
124: PCISSetUp - sets up the `PC_IS` portion of `PCNN` and `PCBDDC` preconditioner context as part of their setup process
126: Input Parameters:
127: + pc - the `PC` object, must be of type `PCNN` or `PCBDDC`
128: . computematrices - Extract the blocks `A_II`, `A_BI`, `A_IB` and `A_BB` from the matrix
129: - computesolvers - Create the `KSP` for the local Dirichlet and Neumann problems
131: Level: advanced
133: .seealso: [](ch_ksp), `PCBDDC`, `PCNN`, `PCISSetUseStiffnessScaling()`, `PCISSetSubdomainDiagonalScaling()`, `PCISScatterArrayNToVecB()`,
134: `PCISSetSubdomainScalingFactor()`,
135: `PCISReset()`, `PCISApplySchur()`, `PCISApplyInvSchur()`
136: @*/
137: PetscErrorCode PCISSetUp(PC pc, PetscBool computematrices, PetscBool computesolvers)
138: {
139: PC_IS *pcis = (PC_IS *)pc->data;
140: Mat_IS *matis;
141: MatReuse reuse;
142: PetscBool flg, issbaij;
144: PetscFunctionBegin;
145: PetscCall(PetscObjectTypeCompare((PetscObject)pc->pmat, MATIS, &flg));
146: PetscCheck(flg, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONG, "Requires matrix for constructing the preconditioner of type MATIS");
147: matis = (Mat_IS *)pc->pmat->data;
148: if (pc->useAmat) {
149: PetscCall(PetscObjectTypeCompare((PetscObject)pc->mat, MATIS, &flg));
150: PetscCheck(flg, 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, *count;
157: PetscInt i;
159: /* get info on mapping */
160: PetscCall(PetscObjectReference((PetscObject)matis->rmapping));
161: PetscCall(ISLocalToGlobalMappingDestroy(&pcis->mapping));
162: pcis->mapping = matis->rmapping;
163: PetscCall(ISLocalToGlobalMappingGetSize(pcis->mapping, &pcis->n));
164: PetscCall(ISLocalToGlobalMappingGetInfo(pcis->mapping, &pcis->n_neigh, &pcis->neigh, &pcis->n_shared, &pcis->shared));
166: /* Identifying interior and interface nodes, in local numbering */
167: PetscCall(ISLocalToGlobalMappingGetNodeInfo(pcis->mapping, NULL, &count, NULL));
168: PetscCall(PetscMalloc1(pcis->n, &idx_I_local));
169: PetscCall(PetscMalloc1(pcis->n, &idx_B_local));
170: for (i = 0, pcis->n_B = 0, n_I = 0; i < pcis->n; i++) {
171: if (count[i] < 2) {
172: idx_I_local[n_I] = i;
173: n_I++;
174: } else {
175: idx_B_local[pcis->n_B] = i;
176: pcis->n_B++;
177: }
178: }
179: PetscCall(ISLocalToGlobalMappingRestoreNodeInfo(pcis->mapping, NULL, &count, NULL));
181: /* Getting the global numbering */
182: idx_B_global = PetscSafePointerPlusOffset(idx_I_local, n_I); /* Just avoiding allocating extra memory, since we have vacant space */
183: idx_I_global = PetscSafePointerPlusOffset(idx_B_local, pcis->n_B);
184: PetscCall(ISLocalToGlobalMappingApply(pcis->mapping, pcis->n_B, idx_B_local, idx_B_global));
185: PetscCall(ISLocalToGlobalMappingApply(pcis->mapping, n_I, idx_I_local, idx_I_global));
187: /* Creating the index sets */
188: PetscCall(ISCreateGeneral(PETSC_COMM_SELF, pcis->n_B, idx_B_local, PETSC_COPY_VALUES, &pcis->is_B_local));
189: PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)pc), pcis->n_B, idx_B_global, PETSC_COPY_VALUES, &pcis->is_B_global));
190: PetscCall(ISCreateGeneral(PETSC_COMM_SELF, n_I, idx_I_local, PETSC_COPY_VALUES, &pcis->is_I_local));
191: PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)pc), n_I, idx_I_global, PETSC_COPY_VALUES, &pcis->is_I_global));
193: /* Freeing memory */
194: PetscCall(PetscFree(idx_B_local));
195: PetscCall(PetscFree(idx_I_local));
197: /* Creating work vectors and arrays */
198: PetscCall(VecDuplicate(matis->x, &pcis->vec1_N));
199: PetscCall(VecDuplicate(pcis->vec1_N, &pcis->vec2_N));
200: PetscCall(VecCreate(PETSC_COMM_SELF, &pcis->vec1_D));
201: PetscCall(VecSetSizes(pcis->vec1_D, pcis->n - pcis->n_B, PETSC_DECIDE));
202: PetscCall(VecSetType(pcis->vec1_D, ((PetscObject)pcis->vec1_N)->type_name));
203: PetscCall(VecDuplicate(pcis->vec1_D, &pcis->vec2_D));
204: PetscCall(VecDuplicate(pcis->vec1_D, &pcis->vec3_D));
205: PetscCall(VecDuplicate(pcis->vec1_D, &pcis->vec4_D));
206: PetscCall(VecCreate(PETSC_COMM_SELF, &pcis->vec1_B));
207: PetscCall(VecSetSizes(pcis->vec1_B, pcis->n_B, PETSC_DECIDE));
208: PetscCall(VecSetType(pcis->vec1_B, ((PetscObject)pcis->vec1_N)->type_name));
209: PetscCall(VecDuplicate(pcis->vec1_B, &pcis->vec2_B));
210: PetscCall(VecDuplicate(pcis->vec1_B, &pcis->vec3_B));
211: PetscCall(MatCreateVecs(pc->pmat, &pcis->vec1_global, NULL));
212: PetscCall(PetscMalloc1(pcis->n, &pcis->work_N));
213: /* scaling vector */
214: if (!pcis->D) { /* it can happen that the user passed in a scaling vector via PCISSetSubdomainDiagonalScaling */
215: PetscCall(VecDuplicate(pcis->vec1_B, &pcis->D));
216: PetscCall(VecSet(pcis->D, pcis->scaling_factor));
217: }
219: /* Creating the scatter contexts */
220: PetscCall(VecScatterCreate(pcis->vec1_N, pcis->is_I_local, pcis->vec1_D, (IS)0, &pcis->N_to_D));
221: PetscCall(VecScatterCreate(pcis->vec1_global, pcis->is_I_global, pcis->vec1_D, (IS)0, &pcis->global_to_D));
222: PetscCall(VecScatterCreate(pcis->vec1_N, pcis->is_B_local, pcis->vec1_B, (IS)0, &pcis->N_to_B));
223: PetscCall(VecScatterCreate(pcis->vec1_global, pcis->is_B_global, pcis->vec1_B, (IS)0, &pcis->global_to_B));
225: /* map from boundary to local */
226: PetscCall(ISLocalToGlobalMappingCreateIS(pcis->is_B_local, &pcis->BtoNmap));
227: }
229: {
230: PetscInt sn;
232: PetscCall(VecGetSize(pcis->D, &sn));
233: if (sn == pcis->n) {
234: PetscCall(VecScatterBegin(pcis->N_to_B, pcis->D, pcis->vec1_B, INSERT_VALUES, SCATTER_FORWARD));
235: PetscCall(VecScatterEnd(pcis->N_to_B, pcis->D, pcis->vec1_B, INSERT_VALUES, SCATTER_FORWARD));
236: PetscCall(VecDestroy(&pcis->D));
237: PetscCall(VecDuplicate(pcis->vec1_B, &pcis->D));
238: PetscCall(VecCopy(pcis->vec1_B, pcis->D));
239: } else PetscCheck(sn == pcis->n_B, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Invalid size for scaling vector. Expected %" PetscInt_FMT " (or full %" PetscInt_FMT "), found %" PetscInt_FMT, pcis->n_B, pcis->n, sn);
240: }
242: /*
243: Extracting the blocks A_II, A_BI, A_IB and A_BB from A. If the numbering
244: is such that interior nodes come first than the interface ones, we have
246: [ A_II | A_IB ]
247: A = [------+------]
248: [ A_BI | A_BB ]
249: */
250: if (computematrices) {
251: PetscBool amat = (PetscBool)(pc->mat != pc->pmat && pc->useAmat);
252: PetscInt bs, ibs;
254: reuse = MAT_INITIAL_MATRIX;
255: if (pcis->reusesubmatrices && pc->setupcalled) {
256: if (pc->flag == SAME_NONZERO_PATTERN) {
257: reuse = MAT_REUSE_MATRIX;
258: } else {
259: reuse = MAT_INITIAL_MATRIX;
260: }
261: }
262: if (reuse == MAT_INITIAL_MATRIX) {
263: PetscCall(MatDestroy(&pcis->A_II));
264: PetscCall(MatDestroy(&pcis->pA_II));
265: PetscCall(MatDestroy(&pcis->A_IB));
266: PetscCall(MatDestroy(&pcis->A_BI));
267: PetscCall(MatDestroy(&pcis->A_BB));
268: }
270: PetscCall(ISLocalToGlobalMappingGetBlockSize(pcis->mapping, &ibs));
271: PetscCall(MatGetBlockSize(matis->A, &bs));
272: PetscCall(MatCreateSubMatrix(matis->A, pcis->is_I_local, pcis->is_I_local, reuse, &pcis->pA_II));
273: if (amat) {
274: Mat_IS *amatis = (Mat_IS *)pc->mat->data;
275: PetscCall(MatCreateSubMatrix(amatis->A, pcis->is_I_local, pcis->is_I_local, reuse, &pcis->A_II));
276: } else {
277: PetscCall(PetscObjectReference((PetscObject)pcis->pA_II));
278: PetscCall(MatDestroy(&pcis->A_II));
279: pcis->A_II = pcis->pA_II;
280: }
281: PetscCall(MatSetBlockSize(pcis->A_II, bs == ibs ? bs : 1));
282: PetscCall(MatSetBlockSize(pcis->pA_II, bs == ibs ? bs : 1));
283: PetscCall(MatCreateSubMatrix(matis->A, pcis->is_B_local, pcis->is_B_local, reuse, &pcis->A_BB));
284: PetscCall(PetscObjectTypeCompare((PetscObject)matis->A, MATSEQSBAIJ, &issbaij));
285: if (!issbaij) {
286: PetscCall(MatCreateSubMatrix(matis->A, pcis->is_I_local, pcis->is_B_local, reuse, &pcis->A_IB));
287: PetscCall(MatCreateSubMatrix(matis->A, pcis->is_B_local, pcis->is_I_local, reuse, &pcis->A_BI));
288: } else {
289: Mat newmat;
291: PetscCall(MatConvert(matis->A, MATSEQBAIJ, MAT_INITIAL_MATRIX, &newmat));
292: PetscCall(MatCreateSubMatrix(newmat, pcis->is_I_local, pcis->is_B_local, reuse, &pcis->A_IB));
293: PetscCall(MatCreateSubMatrix(newmat, pcis->is_B_local, pcis->is_I_local, reuse, &pcis->A_BI));
294: PetscCall(MatDestroy(&newmat));
295: }
296: PetscCall(MatSetBlockSize(pcis->A_BB, bs == ibs ? bs : 1));
297: }
299: /* Creating scaling vector D */
300: PetscCall(PetscOptionsGetBool(((PetscObject)pc)->options, ((PetscObject)pc)->prefix, "-pc_is_use_stiffness_scaling", &pcis->use_stiffness_scaling, NULL));
301: if (pcis->use_stiffness_scaling) {
302: PetscScalar *a;
303: PetscInt i, n;
305: if (pcis->A_BB) {
306: PetscCall(MatGetDiagonal(pcis->A_BB, pcis->D));
307: } else {
308: PetscCall(MatGetDiagonal(matis->A, pcis->vec1_N));
309: PetscCall(VecScatterBegin(pcis->N_to_B, pcis->vec1_N, pcis->D, INSERT_VALUES, SCATTER_FORWARD));
310: PetscCall(VecScatterEnd(pcis->N_to_B, pcis->vec1_N, pcis->D, INSERT_VALUES, SCATTER_FORWARD));
311: }
312: PetscCall(VecAbs(pcis->D));
313: PetscCall(VecGetLocalSize(pcis->D, &n));
314: PetscCall(VecGetArray(pcis->D, &a));
315: for (i = 0; i < n; i++)
316: if (PetscAbsScalar(a[i]) < PETSC_SMALL) a[i] = 1.0;
317: PetscCall(VecRestoreArray(pcis->D, &a));
318: }
319: PetscCall(VecSet(pcis->vec1_global, 0.0));
320: PetscCall(VecScatterBegin(pcis->global_to_B, pcis->D, pcis->vec1_global, ADD_VALUES, SCATTER_REVERSE));
321: PetscCall(VecScatterEnd(pcis->global_to_B, pcis->D, pcis->vec1_global, ADD_VALUES, SCATTER_REVERSE));
322: PetscCall(VecScatterBegin(pcis->global_to_B, pcis->vec1_global, pcis->vec1_B, INSERT_VALUES, SCATTER_FORWARD));
323: PetscCall(VecScatterEnd(pcis->global_to_B, pcis->vec1_global, pcis->vec1_B, INSERT_VALUES, SCATTER_FORWARD));
324: PetscCall(VecPointwiseDivide(pcis->D, pcis->D, pcis->vec1_B));
326: /* Creating the KSP contexts for the local Dirichlet and Neumann problems */
327: if (computesolvers) {
328: PC pc_ctx;
330: pcis->pure_neumann = matis->pure_neumann;
331: /* Dirichlet */
332: PetscCall(KSPCreate(PETSC_COMM_SELF, &pcis->ksp_D));
333: PetscCall(KSPSetNestLevel(pcis->ksp_D, pc->kspnestlevel));
334: PetscCall(KSPSetErrorIfNotConverged(pcis->ksp_D, pc->erroriffailure));
335: PetscCall(PetscObjectIncrementTabLevel((PetscObject)pcis->ksp_D, (PetscObject)pc, 1));
336: PetscCall(KSPSetOperators(pcis->ksp_D, pcis->A_II, pcis->A_II));
337: PetscCall(KSPSetOptionsPrefix(pcis->ksp_D, "is_localD_"));
338: PetscCall(KSPGetPC(pcis->ksp_D, &pc_ctx));
339: PetscCall(PCSetType(pc_ctx, PCLU));
340: PetscCall(KSPSetType(pcis->ksp_D, KSPPREONLY));
341: PetscCall(KSPSetFromOptions(pcis->ksp_D));
342: /* the vectors in the following line are dummy arguments, just telling the KSP the vector size. Values are not used */
343: PetscCall(KSPSetUp(pcis->ksp_D));
344: /* Neumann */
345: PetscCall(KSPCreate(PETSC_COMM_SELF, &pcis->ksp_N));
346: PetscCall(KSPSetNestLevel(pcis->ksp_N, pc->kspnestlevel));
347: PetscCall(KSPSetErrorIfNotConverged(pcis->ksp_N, pc->erroriffailure));
348: PetscCall(PetscObjectIncrementTabLevel((PetscObject)pcis->ksp_N, (PetscObject)pc, 1));
349: PetscCall(KSPSetOperators(pcis->ksp_N, matis->A, matis->A));
350: PetscCall(KSPSetOptionsPrefix(pcis->ksp_N, "is_localN_"));
351: PetscCall(KSPGetPC(pcis->ksp_N, &pc_ctx));
352: PetscCall(PCSetType(pc_ctx, PCLU));
353: PetscCall(KSPSetType(pcis->ksp_N, KSPPREONLY));
354: PetscCall(KSPSetFromOptions(pcis->ksp_N));
355: {
356: 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;
357: PetscReal fixed_factor, floating_factor;
359: PetscCall(PetscOptionsGetReal(((PetscObject)pc_ctx)->options, ((PetscObject)pc_ctx)->prefix, "-pc_is_damp_fixed", &fixed_factor, &damp_fixed));
360: if (!damp_fixed) fixed_factor = 0.0;
361: PetscCall(PetscOptionsGetBool(((PetscObject)pc_ctx)->options, ((PetscObject)pc_ctx)->prefix, "-pc_is_damp_fixed", &damp_fixed, NULL));
363: PetscCall(PetscOptionsGetBool(((PetscObject)pc_ctx)->options, ((PetscObject)pc_ctx)->prefix, "-pc_is_remove_nullspace_fixed", &remove_nullspace_fixed, NULL));
365: PetscCall(PetscOptionsGetReal(((PetscObject)pc_ctx)->options, ((PetscObject)pc_ctx)->prefix, "-pc_is_set_damping_factor_floating", &floating_factor, &set_damping_factor_floating));
366: if (!set_damping_factor_floating) floating_factor = 0.0;
367: PetscCall(PetscOptionsGetBool(((PetscObject)pc_ctx)->options, ((PetscObject)pc_ctx)->prefix, "-pc_is_set_damping_factor_floating", &set_damping_factor_floating, NULL));
368: if (!set_damping_factor_floating) floating_factor = 1.e-12;
370: PetscCall(PetscOptionsGetBool(((PetscObject)pc_ctx)->options, ((PetscObject)pc_ctx)->prefix, "-pc_is_not_damp_floating", ¬_damp_floating, NULL));
372: PetscCall(PetscOptionsGetBool(((PetscObject)pc_ctx)->options, ((PetscObject)pc_ctx)->prefix, "-pc_is_not_remove_nullspace_floating", ¬_remove_nullspace_floating, NULL));
374: if (pcis->pure_neumann) { /* floating subdomain */
375: if (!(not_damp_floating)) {
376: PetscCall(PCFactorSetShiftType(pc_ctx, MAT_SHIFT_NONZERO));
377: PetscCall(PCFactorSetShiftAmount(pc_ctx, floating_factor));
378: }
379: if (!(not_remove_nullspace_floating)) {
380: MatNullSpace nullsp;
381: PetscCall(MatNullSpaceCreate(PETSC_COMM_SELF, PETSC_TRUE, 0, NULL, &nullsp));
382: PetscCall(MatSetNullSpace(matis->A, nullsp));
383: PetscCall(MatNullSpaceDestroy(&nullsp));
384: }
385: } else { /* fixed subdomain */
386: if (damp_fixed) {
387: PetscCall(PCFactorSetShiftType(pc_ctx, MAT_SHIFT_NONZERO));
388: PetscCall(PCFactorSetShiftAmount(pc_ctx, floating_factor));
389: }
390: if (remove_nullspace_fixed) {
391: MatNullSpace nullsp;
392: PetscCall(MatNullSpaceCreate(PETSC_COMM_SELF, PETSC_TRUE, 0, NULL, &nullsp));
393: PetscCall(MatSetNullSpace(matis->A, nullsp));
394: PetscCall(MatNullSpaceDestroy(&nullsp));
395: }
396: }
397: }
398: /* the vectors in the following line are dummy arguments, just telling the KSP the vector size. Values are not used */
399: PetscCall(KSPSetUp(pcis->ksp_N));
400: }
401: PetscFunctionReturn(PETSC_SUCCESS);
402: }
404: /*@
405: PCISReset - Removes all the `PC_IS` parts of the `PC` implementation data structure
407: Input Parameter:
408: . pc - the `PC` object, must be of type `PCNN` or `PCBDDC`
410: Level: advanced
412: .seealso: [](ch_ksp), `PCISSetUseStiffnessScaling()`, `PCISSetSubdomainDiagonalScaling()`, `PCISScatterArrayNToVecB()`, `PCISSetSubdomainScalingFactor()`,
413: `PCISInitialize()`, `PCISApplySchur()`, `PCISApplyInvSchur()`
414: @*/
415: PetscErrorCode PCISReset(PC pc)
416: {
417: PC_IS *pcis = (PC_IS *)pc->data;
418: PetscBool correcttype;
420: PetscFunctionBegin;
421: if (!pc) PetscFunctionReturn(PETSC_SUCCESS);
422: PetscCall(PetscObjectTypeCompareAny((PetscObject)pc, &correcttype, PCBDDC, PCNN, ""));
423: PetscCheck(correcttype, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONGSTATE, "PC must be of type PCNN or PCBDDC");
424: PetscCall(ISDestroy(&pcis->is_B_local));
425: PetscCall(ISDestroy(&pcis->is_I_local));
426: PetscCall(ISDestroy(&pcis->is_B_global));
427: PetscCall(ISDestroy(&pcis->is_I_global));
428: PetscCall(MatDestroy(&pcis->A_II));
429: PetscCall(MatDestroy(&pcis->pA_II));
430: PetscCall(MatDestroy(&pcis->A_IB));
431: PetscCall(MatDestroy(&pcis->A_BI));
432: PetscCall(MatDestroy(&pcis->A_BB));
433: PetscCall(VecDestroy(&pcis->D));
434: PetscCall(KSPDestroy(&pcis->ksp_N));
435: PetscCall(KSPDestroy(&pcis->ksp_D));
436: PetscCall(VecDestroy(&pcis->vec1_N));
437: PetscCall(VecDestroy(&pcis->vec2_N));
438: PetscCall(VecDestroy(&pcis->vec1_D));
439: PetscCall(VecDestroy(&pcis->vec2_D));
440: PetscCall(VecDestroy(&pcis->vec3_D));
441: PetscCall(VecDestroy(&pcis->vec4_D));
442: PetscCall(VecDestroy(&pcis->vec1_B));
443: PetscCall(VecDestroy(&pcis->vec2_B));
444: PetscCall(VecDestroy(&pcis->vec3_B));
445: PetscCall(VecDestroy(&pcis->vec1_global));
446: PetscCall(VecScatterDestroy(&pcis->global_to_D));
447: PetscCall(VecScatterDestroy(&pcis->N_to_B));
448: PetscCall(VecScatterDestroy(&pcis->N_to_D));
449: PetscCall(VecScatterDestroy(&pcis->global_to_B));
450: PetscCall(PetscFree(pcis->work_N));
451: if (pcis->n_neigh > -1) PetscCall(ISLocalToGlobalMappingRestoreInfo(pcis->mapping, &pcis->n_neigh, &pcis->neigh, &pcis->n_shared, &pcis->shared));
452: PetscCall(ISLocalToGlobalMappingDestroy(&pcis->mapping));
453: PetscCall(ISLocalToGlobalMappingDestroy(&pcis->BtoNmap));
454: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCISSetUseStiffnessScaling_C", NULL));
455: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCISSetSubdomainScalingFactor_C", NULL));
456: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCISSetSubdomainDiagonalScaling_C", NULL));
457: PetscFunctionReturn(PETSC_SUCCESS);
458: }
460: /*@
461: PCISInitialize - initializes the `PC_IS` portion of `PCNN` and `PCBDDC` preconditioner context
463: Input Parameter:
464: . pc - the `PC` object, must be of type `PCNN` or `PCBDDC`
466: Level: advanced
468: Note:
469: There is no preconditioner the `PCIS` prefixed routines provide functionality needed by `PCNN` or `PCBDDC`
471: .seealso: [](ch_ksp), `PCBDDC`, `PCNN`, `PCISSetUseStiffnessScaling()`, `PCISSetSubdomainDiagonalScaling()`, `PCISScatterArrayNToVecB()`,
472: `PCISSetSubdomainScalingFactor()`,
473: `PCISReset()`, `PCISApplySchur()`, `PCISApplyInvSchur()`
474: @*/
475: PetscErrorCode PCISInitialize(PC pc)
476: {
477: PC_IS *pcis = (PC_IS *)pc->data;
478: PetscBool correcttype;
480: PetscFunctionBegin;
481: PetscCheck(pcis, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONGSTATE, "PC_IS context must be created by caller");
482: PetscCall(PetscObjectTypeCompareAny((PetscObject)pc, &correcttype, PCBDDC, PCNN, ""));
483: PetscCheck(correcttype, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONGSTATE, "PC must be of type PCNN or PCBDDC");
484: pcis->n_neigh = -1;
485: pcis->scaling_factor = 1.0;
486: pcis->reusesubmatrices = PETSC_TRUE;
487: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCISSetUseStiffnessScaling_C", PCISSetUseStiffnessScaling_IS));
488: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCISSetSubdomainScalingFactor_C", PCISSetSubdomainScalingFactor_IS));
489: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCISSetSubdomainDiagonalScaling_C", PCISSetSubdomainDiagonalScaling_IS));
490: PetscFunctionReturn(PETSC_SUCCESS);
491: }
493: /*@
494: PCISApplySchur - applies the Schur complement arising from the `MATIS` inside the `PCNN` preconditioner
496: Input Parameters:
497: + pc - preconditioner context
498: . v - vector to which the Schur complement is to be applied (it is NOT modified inside this function, UNLESS vec2_B is null)
499: . vec1_B - location to store the result of Schur complement applied to chunk
500: . vec2_B - workspace or `NULL`, `v` is used as workspace in that case
501: . vec1_D - work space
502: - vec2_D - work space
504: Level: advanced
506: .seealso: [](ch_ksp), `PCBDDC`, `PCNN`, `PCISSetUseStiffnessScaling()`, `PCISSetSubdomainDiagonalScaling()`, `PCISScatterArrayNToVecB()`,
507: `PCISSetSubdomainScalingFactor()`, `PCISApplyInvSchur()`,
508: `PCISReset()`, `PCISInitialize()`
509: @*/
510: PetscErrorCode PCISApplySchur(PC pc, Vec v, Vec vec1_B, Vec vec2_B, Vec vec1_D, Vec vec2_D)
511: {
512: PC_IS *pcis = (PC_IS *)pc->data;
514: PetscFunctionBegin;
515: if (!vec2_B) vec2_B = v;
517: PetscCall(MatMult(pcis->A_BB, v, vec1_B));
518: PetscCall(MatMult(pcis->A_IB, v, vec1_D));
519: PetscCall(KSPSolve(pcis->ksp_D, vec1_D, vec2_D));
520: PetscCall(KSPCheckSolve(pcis->ksp_D, pc, vec2_D));
521: PetscCall(MatMult(pcis->A_BI, vec2_D, vec2_B));
522: PetscCall(VecAXPY(vec1_B, -1.0, vec2_B));
523: PetscFunctionReturn(PETSC_SUCCESS);
524: }
526: /*@
527: PCISScatterArrayNToVecB - Scatters interface node values from a big array (of all local nodes, interior or interface,
528: including ghosts) into an interface vector, when in `SCATTER_FORWARD` mode, or vice-versa, when in `SCATTER_REVERSE`
529: mode.
531: Input Parameters:
532: + pc - preconditioner context
533: . array_N - [when in `SCATTER_FORWARD` mode] Array to be scattered into the vector otherwise output array
534: . imode - insert mode, `ADD_VALUES` or `INSERT_VALUES`
535: . smode - scatter mode, `SCATTER_FORWARD` or `SCATTER_REVERSE` mode]
536: - v_B - [when in `SCATTER_REVERSE` mode] Vector to be scattered into the array, otherwise output vector
538: Level: advanced
540: Note:
541: The entries in the array that do not correspond to interface nodes remain unaltered.
543: .seealso: [](ch_ksp), `PCBDDC`, `PCNN`, `PCISSetUseStiffnessScaling()`, `PCISSetSubdomainDiagonalScaling()`,
544: `PCISSetSubdomainScalingFactor()`, `PCISApplySchur()`, `PCISApplyInvSchur()`,
545: `PCISReset()`, `PCISInitialize()`, `InsertMode`
546: @*/
547: PetscErrorCode PCISScatterArrayNToVecB(PC pc, PetscScalar *array_N, Vec v_B, InsertMode imode, ScatterMode smode)
548: {
549: PetscInt i;
550: const PetscInt *idex;
551: PetscScalar *array_B;
552: PC_IS *pcis = (PC_IS *)pc->data;
554: PetscFunctionBegin;
555: PetscCall(VecGetArray(v_B, &array_B));
556: PetscCall(ISGetIndices(pcis->is_B_local, &idex));
558: if (smode == SCATTER_FORWARD) {
559: if (imode == INSERT_VALUES) {
560: for (i = 0; i < pcis->n_B; i++) array_B[i] = array_N[idex[i]];
561: } else { /* ADD_VALUES */
562: for (i = 0; i < pcis->n_B; i++) array_B[i] += array_N[idex[i]];
563: }
564: } else { /* SCATTER_REVERSE */
565: if (imode == INSERT_VALUES) {
566: for (i = 0; i < pcis->n_B; i++) array_N[idex[i]] = array_B[i];
567: } else { /* ADD_VALUES */
568: for (i = 0; i < pcis->n_B; i++) array_N[idex[i]] += array_B[i];
569: }
570: }
571: PetscCall(ISRestoreIndices(pcis->is_B_local, &idex));
572: PetscCall(VecRestoreArray(v_B, &array_B));
573: PetscFunctionReturn(PETSC_SUCCESS);
574: }
576: /*@
577: PCISApplyInvSchur - Solves the Neumann problem related to applying the inverse of the Schur complement.
579: Input Parameters:
580: + pc - preconditioner context
581: . b - vector of local interface nodes (including ghosts)
582: . x - vector of local interface nodes (including ghosts); returns the application of the inverse of the Schur complement to `b`
583: . vec1_N - vector of local nodes (interior and interface, including ghosts); used as work space
584: - vec2_N - vector of local nodes (interior and interface, including ghosts); used as work space
586: Level: advanced
588: Note:
589: Solves the problem
590: .vb
591: [ A_II A_IB ] [ . ] [ 0 ]
592: [ ] [ ] = [ ]
593: [ A_BI A_BB ] [ x ] [ b ]
594: .ve
596: .seealso: [](ch_ksp), `PCBDDC`, `PCNN`, `PCISSetUseStiffnessScaling()`, `PCISSetSubdomainDiagonalScaling()`, `PCISScatterArrayNToVecB()`,
597: `PCISSetSubdomainScalingFactor()`,
598: `PCISReset()`, `PCISInitialize()`
599: @*/
600: PetscErrorCode PCISApplyInvSchur(PC pc, Vec b, Vec x, Vec vec1_N, Vec vec2_N)
601: {
602: PC_IS *pcis = (PC_IS *)pc->data;
604: PetscFunctionBegin;
605: /*
606: Neumann solvers.
607: Applying the inverse of the local Schur complement, i.e, solving a Neumann
608: Problem with zero at the interior nodes of the RHS and extracting the interface
609: part of the solution. inverse Schur complement is applied to b and the result
610: is stored in x.
611: */
612: /* Setting the RHS vec1_N */
613: PetscCall(VecSet(vec1_N, 0.0));
614: PetscCall(VecScatterBegin(pcis->N_to_B, b, vec1_N, INSERT_VALUES, SCATTER_REVERSE));
615: PetscCall(VecScatterEnd(pcis->N_to_B, b, vec1_N, INSERT_VALUES, SCATTER_REVERSE));
616: /* Checking for consistency of the RHS */
617: {
618: PetscBool flg = PETSC_FALSE;
619: PetscCall(PetscOptionsGetBool(NULL, NULL, "-pc_is_check_consistency", &flg, NULL));
620: if (flg) {
621: PetscScalar average;
622: PetscViewer viewer;
623: PetscCall(PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)pc), &viewer));
625: PetscCall(VecSum(vec1_N, &average));
626: average = average / ((PetscReal)pcis->n);
627: PetscCall(PetscViewerASCIIPushSynchronized(viewer));
628: if (pcis->pure_neumann) {
629: PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, "Subdomain %04d is floating. Average = % 1.14e\n", PetscGlobalRank, (double)PetscAbsScalar(average)));
630: } else {
631: PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, "Subdomain %04d is fixed. Average = % 1.14e\n", PetscGlobalRank, (double)PetscAbsScalar(average)));
632: }
633: PetscCall(PetscViewerFlush(viewer));
634: PetscCall(PetscViewerASCIIPopSynchronized(viewer));
635: }
636: }
637: /* Solving the system for vec2_N */
638: PetscCall(KSPSolve(pcis->ksp_N, vec1_N, vec2_N));
639: PetscCall(KSPCheckSolve(pcis->ksp_N, pc, vec2_N));
640: /* Extracting the local interface vector out of the solution */
641: PetscCall(VecScatterBegin(pcis->N_to_B, vec2_N, x, INSERT_VALUES, SCATTER_FORWARD));
642: PetscCall(VecScatterEnd(pcis->N_to_B, vec2_N, x, INSERT_VALUES, SCATTER_FORWARD));
643: PetscFunctionReturn(PETSC_SUCCESS);
644: }