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