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