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