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",&not_damp_floating,NULL);CHKERRQ(ierr);
382 
383       ierr = PetscOptionsGetBool(((PetscObject)pc_ctx)->options,((PetscObject)pc_ctx)->prefix,"-pc_is_not_remove_nullspace_floating",&not_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