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