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