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