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