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