1 #include <petsc/private/pcisimpl.h> /*I "petscpc.h" I*/
2
PCISSetUseStiffnessScaling_IS(PC pc,PetscBool use)3 static PetscErrorCode PCISSetUseStiffnessScaling_IS(PC pc, PetscBool use)
4 {
5 PC_IS *pcis = (PC_IS *)pc->data;
6
7 PetscFunctionBegin;
8 pcis->use_stiffness_scaling = use;
9 PetscFunctionReturn(PETSC_SUCCESS);
10 }
11
12 /*@
13 PCISSetUseStiffnessScaling - Tells `PCIS` to construct partition of unity using
14 the local matrices' diagonal entries
15
16 Logically Collective
17
18 Input Parameters:
19 + pc - the preconditioning context
20 - use - whether or not it should use matrix diagonal to build partition of unity.
21
22 Level: intermediate
23
24 .seealso: [](ch_ksp), `PCBDDC`, `PCNN`, `PCISSetSubdomainDiagonalScaling()`, `PCISScatterArrayNToVecB()`,
25 `PCISSetSubdomainScalingFactor()`,
26 `PCISReset()`, `PCISInitialize()`, `PCISApplyInvSchur()`, `PCISApplySchur()`
27 @*/
PCISSetUseStiffnessScaling(PC pc,PetscBool use)28 PetscErrorCode PCISSetUseStiffnessScaling(PC pc, PetscBool use)
29 {
30 PetscFunctionBegin;
31 PetscValidHeaderSpecific(pc, PC_CLASSID, 1);
32 PetscValidLogicalCollectiveBool(pc, use, 2);
33 PetscTryMethod(pc, "PCISSetUseStiffnessScaling_C", (PC, PetscBool), (pc, use));
34 PetscFunctionReturn(PETSC_SUCCESS);
35 }
36
PCISSetSubdomainDiagonalScaling_IS(PC pc,Vec scaling_factors)37 static PetscErrorCode PCISSetSubdomainDiagonalScaling_IS(PC pc, Vec scaling_factors)
38 {
39 PC_IS *pcis = (PC_IS *)pc->data;
40
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;
47
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 }
59
60 /*@
61 PCISSetSubdomainDiagonalScaling - Set diagonal scaling for `PCIS`.
62
63 Logically Collective
64
65 Input Parameters:
66 + pc - the preconditioning context
67 - scaling_factors - scaling factors for the subdomain
68
69 Level: intermediate
70
71 Note:
72 Intended for use with jumping coefficients cases.
73
74 .seealso: [](ch_ksp), `PCBDDC`, `PCNN`, `PCISScatterArrayNToVecB()`,
75 `PCISSetSubdomainScalingFactor()`, `PCISSetUseStiffnessScaling()`,
76 `PCISReset()`, `PCISInitialize()`, `PCISApplyInvSchur()`, `PCISApplySchur()`
77 @*/
PCISSetSubdomainDiagonalScaling(PC pc,Vec scaling_factors)78 PetscErrorCode PCISSetSubdomainDiagonalScaling(PC pc, Vec scaling_factors)
79 {
80 PetscFunctionBegin;
81 PetscValidHeaderSpecific(pc, PC_CLASSID, 1);
82 PetscValidHeaderSpecific(scaling_factors, VEC_CLASSID, 2);
83 PetscTryMethod(pc, "PCISSetSubdomainDiagonalScaling_C", (PC, Vec), (pc, scaling_factors));
84 PetscFunctionReturn(PETSC_SUCCESS);
85 }
86
PCISSetSubdomainScalingFactor_IS(PC pc,PetscScalar scal)87 static PetscErrorCode PCISSetSubdomainScalingFactor_IS(PC pc, PetscScalar scal)
88 {
89 PC_IS *pcis = (PC_IS *)pc->data;
90
91 PetscFunctionBegin;
92 pcis->scaling_factor = scal;
93 if (pcis->D) PetscCall(VecSet(pcis->D, pcis->scaling_factor));
94 PetscFunctionReturn(PETSC_SUCCESS);
95 }
96
97 /*@
98 PCISSetSubdomainScalingFactor - Set scaling factor for `PCIS`.
99
100 Not Collective
101
102 Input Parameters:
103 + pc - the preconditioning context
104 - scal - scaling factor for the subdomain
105
106 Level: intermediate
107
108 Note:
109 Intended for use with the jumping coefficients cases.
110
111 .seealso: [](ch_ksp), `PCBDDC`, `PCNN`, `PCISScatterArrayNToVecB()`,
112 `PCISSetSubdomainDiagonalScaling()`, `PCISSetUseStiffnessScaling()`,
113 `PCISReset()`, `PCISInitialize()`, `PCISApplyInvSchur()`, `PCISApplySchur()`
114 @*/
PCISSetSubdomainScalingFactor(PC pc,PetscScalar scal)115 PetscErrorCode PCISSetSubdomainScalingFactor(PC pc, PetscScalar scal)
116 {
117 PetscFunctionBegin;
118 PetscValidHeaderSpecific(pc, PC_CLASSID, 1);
119 PetscTryMethod(pc, "PCISSetSubdomainScalingFactor_C", (PC, PetscScalar), (pc, scal));
120 PetscFunctionReturn(PETSC_SUCCESS);
121 }
122
123 /*@
124 PCISSetUp - sets up the `PC_IS` portion of `PCNN` and `PCBDDC` preconditioner context as part of their setup process
125
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
130
131 Level: advanced
132
133 .seealso: [](ch_ksp), `PCBDDC`, `PCNN`, `PCISSetUseStiffnessScaling()`, `PCISSetSubdomainDiagonalScaling()`, `PCISScatterArrayNToVecB()`,
134 `PCISSetSubdomainScalingFactor()`,
135 `PCISReset()`, `PCISApplySchur()`, `PCISApplyInvSchur()`
136 @*/
PCISSetUp(PC pc,PetscBool computematrices,PetscBool computesolvers)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;
143
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 }
152
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;
158
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));
165
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));
180
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));
186
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));
192
193 /* Freeing memory */
194 PetscCall(PetscFree(idx_B_local));
195 PetscCall(PetscFree(idx_I_local));
196
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 }
218
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));
224
225 /* map from boundary to local */
226 PetscCall(ISLocalToGlobalMappingCreateIS(pcis->is_B_local, &pcis->BtoNmap));
227 }
228
229 {
230 PetscInt sn;
231
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 }
241
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
245
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;
253
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 }
269
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;
290
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 }
298
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;
304
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
326 /* Creating the KSP contexts for the local Dirichlet and Neumann problems */
327 if (computesolvers) {
328 PC pc_ctx;
329
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;
358
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));
362
363 PetscCall(PetscOptionsGetBool(((PetscObject)pc_ctx)->options, ((PetscObject)pc_ctx)->prefix, "-pc_is_remove_nullspace_fixed", &remove_nullspace_fixed, NULL));
364
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;
369
370 PetscCall(PetscOptionsGetBool(((PetscObject)pc_ctx)->options, ((PetscObject)pc_ctx)->prefix, "-pc_is_not_damp_floating", ¬_damp_floating, NULL));
371
372 PetscCall(PetscOptionsGetBool(((PetscObject)pc_ctx)->options, ((PetscObject)pc_ctx)->prefix, "-pc_is_not_remove_nullspace_floating", ¬_remove_nullspace_floating, NULL));
373
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 }
403
404 /*@
405 PCISReset - Removes all the `PC_IS` parts of the `PC` implementation data structure
406
407 Input Parameter:
408 . pc - the `PC` object, must be of type `PCNN` or `PCBDDC`
409
410 Level: advanced
411
412 .seealso: [](ch_ksp), `PCISSetUseStiffnessScaling()`, `PCISSetSubdomainDiagonalScaling()`, `PCISScatterArrayNToVecB()`, `PCISSetSubdomainScalingFactor()`,
413 `PCISInitialize()`, `PCISApplySchur()`, `PCISApplyInvSchur()`
414 @*/
PCISReset(PC pc)415 PetscErrorCode PCISReset(PC pc)
416 {
417 PC_IS *pcis = (PC_IS *)pc->data;
418 PetscBool correcttype;
419
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 }
459
460 /*@
461 PCISInitialize - initializes the `PC_IS` portion of `PCNN` and `PCBDDC` preconditioner context
462
463 Input Parameter:
464 . pc - the `PC` object, must be of type `PCNN` or `PCBDDC`
465
466 Level: advanced
467
468 Note:
469 There is no preconditioner the `PCIS` prefixed routines provide functionality needed by `PCNN` or `PCBDDC`
470
471 .seealso: [](ch_ksp), `PCBDDC`, `PCNN`, `PCISSetUseStiffnessScaling()`, `PCISSetSubdomainDiagonalScaling()`, `PCISScatterArrayNToVecB()`,
472 `PCISSetSubdomainScalingFactor()`,
473 `PCISReset()`, `PCISApplySchur()`, `PCISApplyInvSchur()`
474 @*/
PCISInitialize(PC pc)475 PetscErrorCode PCISInitialize(PC pc)
476 {
477 PC_IS *pcis = (PC_IS *)pc->data;
478 PetscBool correcttype;
479
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 }
492
493 /*@
494 PCISApplySchur - applies the Schur complement arising from the `MATIS` inside the `PCNN` preconditioner
495
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
503
504 Level: advanced
505
506 .seealso: [](ch_ksp), `PCBDDC`, `PCNN`, `PCISSetUseStiffnessScaling()`, `PCISSetSubdomainDiagonalScaling()`, `PCISScatterArrayNToVecB()`,
507 `PCISSetSubdomainScalingFactor()`, `PCISApplyInvSchur()`,
508 `PCISReset()`, `PCISInitialize()`
509 @*/
PCISApplySchur(PC pc,Vec v,Vec vec1_B,Vec vec2_B,Vec vec1_D,Vec vec2_D)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;
513
514 PetscFunctionBegin;
515 if (!vec2_B) vec2_B = v;
516
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 }
525
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.
530
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
537
538 Level: advanced
539
540 Note:
541 The entries in the array that do not correspond to interface nodes remain unaltered.
542
543 .seealso: [](ch_ksp), `PCBDDC`, `PCNN`, `PCISSetUseStiffnessScaling()`, `PCISSetSubdomainDiagonalScaling()`,
544 `PCISSetSubdomainScalingFactor()`, `PCISApplySchur()`, `PCISApplyInvSchur()`,
545 `PCISReset()`, `PCISInitialize()`, `InsertMode`
546 @*/
PCISScatterArrayNToVecB(PC pc,PetscScalar * array_N,Vec v_B,InsertMode imode,ScatterMode smode)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;
553
554 PetscFunctionBegin;
555 PetscCall(VecGetArray(v_B, &array_B));
556 PetscCall(ISGetIndices(pcis->is_B_local, &idex));
557
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 }
575
576 /*@
577 PCISApplyInvSchur - Solves the Neumann problem related to applying the inverse of the Schur complement.
578
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
585
586 Level: advanced
587
588 Note:
589 Solves the problem
590 .vb
591 [ A_II A_IB ] [ . ] [ 0 ]
592 [ ] [ ] = [ ]
593 [ A_BI A_BB ] [ x ] [ b ]
594 .ve
595
596 .seealso: [](ch_ksp), `PCBDDC`, `PCNN`, `PCISSetUseStiffnessScaling()`, `PCISSetSubdomainDiagonalScaling()`, `PCISScatterArrayNToVecB()`,
597 `PCISSetSubdomainScalingFactor()`,
598 `PCISReset()`, `PCISInitialize()`
599 @*/
PCISApplyInvSchur(PC pc,Vec b,Vec x,Vec vec1_N,Vec vec2_N)600 PetscErrorCode PCISApplyInvSchur(PC pc, Vec b, Vec x, Vec vec1_N, Vec vec2_N)
601 {
602 PC_IS *pcis = (PC_IS *)pc->data;
603
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));
624
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 }
645