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 preconditioning matrix 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));
325: /* See historical note 01, at the bottom of this file. */
327: /* Creating the KSP contexts for the local Dirichlet and Neumann problems */
328: if (computesolvers) {
329: PC pc_ctx;
331: pcis->pure_neumann = matis->pure_neumann;
332: /* Dirichlet */
333: PetscCall(KSPCreate(PETSC_COMM_SELF, &pcis->ksp_D));
334: PetscCall(KSPSetNestLevel(pcis->ksp_D, pc->kspnestlevel));
335: PetscCall(KSPSetErrorIfNotConverged(pcis->ksp_D, pc->erroriffailure));
336: PetscCall(PetscObjectIncrementTabLevel((PetscObject)pcis->ksp_D, (PetscObject)pc, 1));
337: PetscCall(KSPSetOperators(pcis->ksp_D, pcis->A_II, pcis->A_II));
338: PetscCall(KSPSetOptionsPrefix(pcis->ksp_D, "is_localD_"));
339: PetscCall(KSPGetPC(pcis->ksp_D, &pc_ctx));
340: PetscCall(PCSetType(pc_ctx, PCLU));
341: PetscCall(KSPSetType(pcis->ksp_D, KSPPREONLY));
342: PetscCall(KSPSetFromOptions(pcis->ksp_D));
343: /* the vectors in the following line are dummy arguments, just telling the KSP the vector size. Values are not used */
344: PetscCall(KSPSetUp(pcis->ksp_D));
345: /* Neumann */
346: PetscCall(KSPCreate(PETSC_COMM_SELF, &pcis->ksp_N));
347: PetscCall(KSPSetNestLevel(pcis->ksp_N, pc->kspnestlevel));
348: PetscCall(KSPSetErrorIfNotConverged(pcis->ksp_N, pc->erroriffailure));
349: PetscCall(PetscObjectIncrementTabLevel((PetscObject)pcis->ksp_N, (PetscObject)pc, 1));
350: PetscCall(KSPSetOperators(pcis->ksp_N, matis->A, matis->A));
351: PetscCall(KSPSetOptionsPrefix(pcis->ksp_N, "is_localN_"));
352: PetscCall(KSPGetPC(pcis->ksp_N, &pc_ctx));
353: PetscCall(PCSetType(pc_ctx, PCLU));
354: PetscCall(KSPSetType(pcis->ksp_N, KSPPREONLY));
355: PetscCall(KSPSetFromOptions(pcis->ksp_N));
356: {
357: 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;
358: PetscReal fixed_factor, floating_factor;
360: PetscCall(PetscOptionsGetReal(((PetscObject)pc_ctx)->options, ((PetscObject)pc_ctx)->prefix, "-pc_is_damp_fixed", &fixed_factor, &damp_fixed));
361: if (!damp_fixed) fixed_factor = 0.0;
362: PetscCall(PetscOptionsGetBool(((PetscObject)pc_ctx)->options, ((PetscObject)pc_ctx)->prefix, "-pc_is_damp_fixed", &damp_fixed, NULL));
364: PetscCall(PetscOptionsGetBool(((PetscObject)pc_ctx)->options, ((PetscObject)pc_ctx)->prefix, "-pc_is_remove_nullspace_fixed", &remove_nullspace_fixed, NULL));
366: PetscCall(PetscOptionsGetReal(((PetscObject)pc_ctx)->options, ((PetscObject)pc_ctx)->prefix, "-pc_is_set_damping_factor_floating", &floating_factor, &set_damping_factor_floating));
367: if (!set_damping_factor_floating) floating_factor = 0.0;
368: PetscCall(PetscOptionsGetBool(((PetscObject)pc_ctx)->options, ((PetscObject)pc_ctx)->prefix, "-pc_is_set_damping_factor_floating", &set_damping_factor_floating, NULL));
369: if (!set_damping_factor_floating) floating_factor = 1.e-12;
371: PetscCall(PetscOptionsGetBool(((PetscObject)pc_ctx)->options, ((PetscObject)pc_ctx)->prefix, "-pc_is_not_damp_floating", ¬_damp_floating, NULL));
373: PetscCall(PetscOptionsGetBool(((PetscObject)pc_ctx)->options, ((PetscObject)pc_ctx)->prefix, "-pc_is_not_remove_nullspace_floating", ¬_remove_nullspace_floating, NULL));
375: if (pcis->pure_neumann) { /* floating subdomain */
376: if (!(not_damp_floating)) {
377: PetscCall(PCFactorSetShiftType(pc_ctx, MAT_SHIFT_NONZERO));
378: PetscCall(PCFactorSetShiftAmount(pc_ctx, floating_factor));
379: }
380: if (!(not_remove_nullspace_floating)) {
381: MatNullSpace nullsp;
382: PetscCall(MatNullSpaceCreate(PETSC_COMM_SELF, PETSC_TRUE, 0, NULL, &nullsp));
383: PetscCall(MatSetNullSpace(matis->A, nullsp));
384: PetscCall(MatNullSpaceDestroy(&nullsp));
385: }
386: } else { /* fixed subdomain */
387: if (damp_fixed) {
388: PetscCall(PCFactorSetShiftType(pc_ctx, MAT_SHIFT_NONZERO));
389: PetscCall(PCFactorSetShiftAmount(pc_ctx, floating_factor));
390: }
391: if (remove_nullspace_fixed) {
392: MatNullSpace nullsp;
393: PetscCall(MatNullSpaceCreate(PETSC_COMM_SELF, PETSC_TRUE, 0, NULL, &nullsp));
394: PetscCall(MatSetNullSpace(matis->A, nullsp));
395: PetscCall(MatNullSpaceDestroy(&nullsp));
396: }
397: }
398: }
399: /* the vectors in the following line are dummy arguments, just telling the KSP the vector size. Values are not used */
400: PetscCall(KSPSetUp(pcis->ksp_N));
401: }
402: PetscFunctionReturn(PETSC_SUCCESS);
403: }
405: /*@
406: PCISReset - Removes all the `PC_IS` parts of the `PC` implementation data structure
408: Input Parameter:
409: . pc - the `PC` object, must be of type `PCNN` or `PCBDDC`
411: Level: advanced
413: .seealso: [](ch_ksp), `PCISSetUseStiffnessScaling()`, `PCISSetSubdomainDiagonalScaling()`, `PCISScatterArrayNToVecB()`, `PCISSetSubdomainScalingFactor()`,
414: `PCISInitialize()`, `PCISApplySchur()`, `PCISApplyInvSchur()`
415: @*/
416: PetscErrorCode PCISReset(PC pc)
417: {
418: PC_IS *pcis = (PC_IS *)pc->data;
419: PetscBool correcttype;
421: PetscFunctionBegin;
422: if (!pc) PetscFunctionReturn(PETSC_SUCCESS);
423: PetscCall(PetscObjectTypeCompareAny((PetscObject)pc, &correcttype, PCBDDC, PCNN, ""));
424: PetscCheck(correcttype, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONGSTATE, "PC must be of type PCNN or PCBDDC");
425: PetscCall(ISDestroy(&pcis->is_B_local));
426: PetscCall(ISDestroy(&pcis->is_I_local));
427: PetscCall(ISDestroy(&pcis->is_B_global));
428: PetscCall(ISDestroy(&pcis->is_I_global));
429: PetscCall(MatDestroy(&pcis->A_II));
430: PetscCall(MatDestroy(&pcis->pA_II));
431: PetscCall(MatDestroy(&pcis->A_IB));
432: PetscCall(MatDestroy(&pcis->A_BI));
433: PetscCall(MatDestroy(&pcis->A_BB));
434: PetscCall(VecDestroy(&pcis->D));
435: PetscCall(KSPDestroy(&pcis->ksp_N));
436: PetscCall(KSPDestroy(&pcis->ksp_D));
437: PetscCall(VecDestroy(&pcis->vec1_N));
438: PetscCall(VecDestroy(&pcis->vec2_N));
439: PetscCall(VecDestroy(&pcis->vec1_D));
440: PetscCall(VecDestroy(&pcis->vec2_D));
441: PetscCall(VecDestroy(&pcis->vec3_D));
442: PetscCall(VecDestroy(&pcis->vec4_D));
443: PetscCall(VecDestroy(&pcis->vec1_B));
444: PetscCall(VecDestroy(&pcis->vec2_B));
445: PetscCall(VecDestroy(&pcis->vec3_B));
446: PetscCall(VecDestroy(&pcis->vec1_global));
447: PetscCall(VecScatterDestroy(&pcis->global_to_D));
448: PetscCall(VecScatterDestroy(&pcis->N_to_B));
449: PetscCall(VecScatterDestroy(&pcis->N_to_D));
450: PetscCall(VecScatterDestroy(&pcis->global_to_B));
451: PetscCall(PetscFree(pcis->work_N));
452: if (pcis->n_neigh > -1) PetscCall(ISLocalToGlobalMappingRestoreInfo(pcis->mapping, &pcis->n_neigh, &pcis->neigh, &pcis->n_shared, &pcis->shared));
453: PetscCall(ISLocalToGlobalMappingDestroy(&pcis->mapping));
454: PetscCall(ISLocalToGlobalMappingDestroy(&pcis->BtoNmap));
455: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCISSetUseStiffnessScaling_C", NULL));
456: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCISSetSubdomainScalingFactor_C", NULL));
457: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCISSetSubdomainDiagonalScaling_C", NULL));
458: PetscFunctionReturn(PETSC_SUCCESS);
459: }
461: /*@
462: PCISInitialize - initializes the `PC_IS` portion of `PCNN` and `PCBDDC` preconditioner context
464: Input Parameter:
465: . pc - the `PC` object, must be of type `PCNN` or `PCBDDC`
467: Level: advanced
469: Note:
470: There is no preconditioner the `PCIS` prefixed routines provide functionality needed by `PCNN` or `PCBDDC`
472: .seealso: [](ch_ksp), `PCBDDC`, `PCNN`, `PCISSetUseStiffnessScaling()`, `PCISSetSubdomainDiagonalScaling()`, `PCISScatterArrayNToVecB()`,
473: `PCISSetSubdomainScalingFactor()`,
474: `PCISReset()`, `PCISApplySchur()`, `PCISApplyInvSchur()`
475: @*/
476: PetscErrorCode PCISInitialize(PC pc)
477: {
478: PC_IS *pcis = (PC_IS *)pc->data;
479: PetscBool correcttype;
481: PetscFunctionBegin;
482: PetscCheck(pcis, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONGSTATE, "PC_IS context must be created by caller");
483: PetscCall(PetscObjectTypeCompareAny((PetscObject)pc, &correcttype, PCBDDC, PCNN, ""));
484: PetscCheck(correcttype, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONGSTATE, "PC must be of type PCNN or PCBDDC");
485: pcis->n_neigh = -1;
486: pcis->scaling_factor = 1.0;
487: pcis->reusesubmatrices = PETSC_TRUE;
488: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCISSetUseStiffnessScaling_C", PCISSetUseStiffnessScaling_IS));
489: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCISSetSubdomainScalingFactor_C", PCISSetSubdomainScalingFactor_IS));
490: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCISSetSubdomainDiagonalScaling_C", PCISSetSubdomainDiagonalScaling_IS));
491: PetscFunctionReturn(PETSC_SUCCESS);
492: }
494: /*@
495: PCISApplySchur - applies the Schur complement arising from the `MATIS` inside the `PCNN` preconditioner
497: Input Parameters:
498: + pc - preconditioner context
499: . v - vector to which the Schur complement is to be applied (it is NOT modified inside this function, UNLESS vec2_B is null)
500: . vec1_B - location to store the result of Schur complement applied to chunk
501: . vec2_B - workspace or `NULL`, `v` is used as workspace in that case
502: . vec1_D - work space
503: - vec2_D - work space
505: Level: advanced
507: .seealso: [](ch_ksp), `PCBDDC`, `PCNN`, `PCISSetUseStiffnessScaling()`, `PCISSetSubdomainDiagonalScaling()`, `PCISScatterArrayNToVecB()`,
508: `PCISSetSubdomainScalingFactor()`, `PCISApplyInvSchur()`,
509: `PCISReset()`, `PCISInitialize()`
510: @*/
511: PetscErrorCode PCISApplySchur(PC pc, Vec v, Vec vec1_B, Vec vec2_B, Vec vec1_D, Vec vec2_D)
512: {
513: PC_IS *pcis = (PC_IS *)pc->data;
515: PetscFunctionBegin;
516: if (!vec2_B) vec2_B = v;
518: PetscCall(MatMult(pcis->A_BB, v, vec1_B));
519: PetscCall(MatMult(pcis->A_IB, v, vec1_D));
520: PetscCall(KSPSolve(pcis->ksp_D, vec1_D, vec2_D));
521: PetscCall(KSPCheckSolve(pcis->ksp_D, pc, vec2_D));
522: PetscCall(MatMult(pcis->A_BI, vec2_D, vec2_B));
523: PetscCall(VecAXPY(vec1_B, -1.0, vec2_B));
524: PetscFunctionReturn(PETSC_SUCCESS);
525: }
527: /*@
528: PCISScatterArrayNToVecB - Scatters interface node values from a big array (of all local nodes, interior or interface,
529: including ghosts) into an interface vector, when in `SCATTER_FORWARD` mode, or vice-versa, when in `SCATTER_REVERSE`
530: mode.
532: Input Parameters:
533: + pc - preconditioner context
534: . array_N - [when in `SCATTER_FORWARD` mode] Array to be scattered into the vector otherwise output array
535: . imode - insert mode, `ADD_VALUES` or `INSERT_VALUES`
536: . smode - scatter mode, `SCATTER_FORWARD` or `SCATTER_REVERSE` mode]
537: - v_B - [when in `SCATTER_REVERSE` mode] Vector to be scattered into the array, otherwise output vector
539: Level: advanced
541: Note:
542: The entries in the array that do not correspond to interface nodes remain unaltered.
544: .seealso: [](ch_ksp), `PCBDDC`, `PCNN`, `PCISSetUseStiffnessScaling()`, `PCISSetSubdomainDiagonalScaling()`,
545: `PCISSetSubdomainScalingFactor()`, `PCISApplySchur()`, `PCISApplyInvSchur()`,
546: `PCISReset()`, `PCISInitialize()`, `InsertMode`
547: @*/
548: PetscErrorCode PCISScatterArrayNToVecB(PC pc, PetscScalar *array_N, Vec v_B, InsertMode imode, ScatterMode smode)
549: {
550: PetscInt i;
551: const PetscInt *idex;
552: PetscScalar *array_B;
553: PC_IS *pcis = (PC_IS *)pc->data;
555: PetscFunctionBegin;
556: PetscCall(VecGetArray(v_B, &array_B));
557: PetscCall(ISGetIndices(pcis->is_B_local, &idex));
559: if (smode == SCATTER_FORWARD) {
560: if (imode == INSERT_VALUES) {
561: for (i = 0; i < pcis->n_B; i++) array_B[i] = array_N[idex[i]];
562: } else { /* ADD_VALUES */
563: for (i = 0; i < pcis->n_B; i++) array_B[i] += array_N[idex[i]];
564: }
565: } else { /* SCATTER_REVERSE */
566: if (imode == INSERT_VALUES) {
567: for (i = 0; i < pcis->n_B; i++) array_N[idex[i]] = array_B[i];
568: } else { /* ADD_VALUES */
569: for (i = 0; i < pcis->n_B; i++) array_N[idex[i]] += array_B[i];
570: }
571: }
572: PetscCall(ISRestoreIndices(pcis->is_B_local, &idex));
573: PetscCall(VecRestoreArray(v_B, &array_B));
574: PetscFunctionReturn(PETSC_SUCCESS);
575: }
577: /*@
578: PCISApplyInvSchur - Solves the Neumann problem related to applying the inverse of the Schur complement.
580: Input Parameters:
581: + pc - preconditioner context
582: . b - vector of local interface nodes (including ghosts)
583: . x - vector of local interface nodes (including ghosts); returns the application of the inverse of the Schur complement to `b`
584: . vec1_N - vector of local nodes (interior and interface, including ghosts); used as work space
585: - vec2_N - vector of local nodes (interior and interface, including ghosts); used as work space
587: Level: advanced
589: Note:
590: Solves the problem
591: .vb
592: [ A_II A_IB ] [ . ] [ 0 ]
593: [ ] [ ] = [ ]
594: [ A_BI A_BB ] [ x ] [ b ]
595: .ve
597: .seealso: [](ch_ksp), `PCBDDC`, `PCNN`, `PCISSetUseStiffnessScaling()`, `PCISSetSubdomainDiagonalScaling()`, `PCISScatterArrayNToVecB()`,
598: `PCISSetSubdomainScalingFactor()`,
599: `PCISReset()`, `PCISInitialize()`
600: @*/
601: PetscErrorCode PCISApplyInvSchur(PC pc, Vec b, Vec x, Vec vec1_N, Vec vec2_N)
602: {
603: PC_IS *pcis = (PC_IS *)pc->data;
605: PetscFunctionBegin;
606: /*
607: Neumann solvers.
608: Applying the inverse of the local Schur complement, i.e, solving a Neumann
609: Problem with zero at the interior nodes of the RHS and extracting the interface
610: part of the solution. inverse Schur complement is applied to b and the result
611: is stored in x.
612: */
613: /* Setting the RHS vec1_N */
614: PetscCall(VecSet(vec1_N, 0.0));
615: PetscCall(VecScatterBegin(pcis->N_to_B, b, vec1_N, INSERT_VALUES, SCATTER_REVERSE));
616: PetscCall(VecScatterEnd(pcis->N_to_B, b, vec1_N, INSERT_VALUES, SCATTER_REVERSE));
617: /* Checking for consistency of the RHS */
618: {
619: PetscBool flg = PETSC_FALSE;
620: PetscCall(PetscOptionsGetBool(NULL, NULL, "-pc_is_check_consistency", &flg, NULL));
621: if (flg) {
622: PetscScalar average;
623: PetscViewer viewer;
624: PetscCall(PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)pc), &viewer));
626: PetscCall(VecSum(vec1_N, &average));
627: average = average / ((PetscReal)pcis->n);
628: PetscCall(PetscViewerASCIIPushSynchronized(viewer));
629: if (pcis->pure_neumann) {
630: PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, "Subdomain %04d is floating. Average = % 1.14e\n", PetscGlobalRank, (double)PetscAbsScalar(average)));
631: } else {
632: PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, "Subdomain %04d is fixed. Average = % 1.14e\n", PetscGlobalRank, (double)PetscAbsScalar(average)));
633: }
634: PetscCall(PetscViewerFlush(viewer));
635: PetscCall(PetscViewerASCIIPopSynchronized(viewer));
636: }
637: }
638: /* Solving the system for vec2_N */
639: PetscCall(KSPSolve(pcis->ksp_N, vec1_N, vec2_N));
640: PetscCall(KSPCheckSolve(pcis->ksp_N, pc, vec2_N));
641: /* Extracting the local interface vector out of the solution */
642: PetscCall(VecScatterBegin(pcis->N_to_B, vec2_N, x, INSERT_VALUES, SCATTER_FORWARD));
643: PetscCall(VecScatterEnd(pcis->N_to_B, vec2_N, x, INSERT_VALUES, SCATTER_FORWARD));
644: PetscFunctionReturn(PETSC_SUCCESS);
645: }