1 /*$Id: is.c,v 1.9 2001/08/07 03:03:41 balay Exp $*/ 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 = VecDuplicate(pc->vec,&counter);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 { 106 Vec global; 107 ierr = PCGetVector(pc,&global);CHKERRQ(ierr); 108 ierr = VecDuplicate(global,&pcis->vec1_global);CHKERRQ(ierr); 109 } 110 ierr = PetscMalloc((pcis->n)*sizeof(PetscScalar),&pcis->work_N);CHKERRQ(ierr); 111 112 /* Creating the scatter contexts */ 113 ierr = VecScatterCreate(pc->vec,pcis->is_I_global,pcis->vec1_D,(IS)0,&pcis->global_to_D);CHKERRQ(ierr); 114 ierr = VecScatterCreate(pcis->vec1_N,pcis->is_B_local,pcis->vec1_B,(IS)0,&pcis->N_to_B);CHKERRQ(ierr); 115 ierr = VecScatterCreate(pc->vec,pcis->is_B_global,pcis->vec1_B,(IS)0,&pcis->global_to_B);CHKERRQ(ierr); 116 117 /* Creating scaling "matrix" D, from information in vec1_N */ 118 ierr = VecDuplicate(pcis->vec1_B,&pcis->D);CHKERRQ(ierr); 119 ierr = VecScatterBegin(pcis->vec1_N,pcis->D,INSERT_VALUES,SCATTER_FORWARD,pcis->N_to_B);CHKERRQ(ierr); 120 ierr = VecScatterEnd (pcis->vec1_N,pcis->D,INSERT_VALUES,SCATTER_FORWARD,pcis->N_to_B);CHKERRQ(ierr); 121 ierr = VecReciprocal(pcis->D);CHKERRQ(ierr); 122 123 /* See historical note 01, at the bottom of this file. */ 124 125 /* 126 Creating the KSP contexts for the local Dirichlet and Neumann problems. 127 */ 128 { 129 PC pc_ctx; 130 /* Dirichlet */ 131 ierr = KSPCreate(PETSC_COMM_SELF,&pcis->ksp_D);CHKERRQ(ierr); 132 ierr = KSPSetOperators(pcis->ksp_D,pcis->A_II,pcis->A_II,SAME_PRECONDITIONER);CHKERRQ(ierr); 133 ierr = KSPSetOptionsPrefix(pcis->ksp_D,"is_localD_");CHKERRQ(ierr); 134 ierr = KSPGetPC(pcis->ksp_D,&pc_ctx);CHKERRQ(ierr); 135 ierr = PCSetType(pc_ctx,PCLU);CHKERRQ(ierr); 136 ierr = KSPSetType(pcis->ksp_D,KSPPREONLY);CHKERRQ(ierr); 137 ierr = KSPSetFromOptions(pcis->ksp_D);CHKERRQ(ierr); 138 /* the vectors in the following line are dummy arguments, just telling the KSP the vector size. Values are not used */ 139 ierr = KSPSetRhs(pcis->ksp_D,pcis->vec1_D);CHKERRQ(ierr); 140 ierr = KSPSetSolution(pcis->ksp_D,pcis->vec2_D);CHKERRQ(ierr); 141 ierr = KSPSetUp(pcis->ksp_D);CHKERRQ(ierr); 142 /* Neumann */ 143 ierr = KSPCreate(PETSC_COMM_SELF,&pcis->ksp_N);CHKERRQ(ierr); 144 ierr = KSPSetOperators(pcis->ksp_N,matis->A,matis->A,SAME_PRECONDITIONER);CHKERRQ(ierr); 145 ierr = KSPSetOptionsPrefix(pcis->ksp_N,"is_localN_");CHKERRQ(ierr); 146 ierr = KSPGetPC(pcis->ksp_N,&pc_ctx);CHKERRQ(ierr); 147 ierr = PCSetType(pc_ctx,PCLU);CHKERRQ(ierr); 148 ierr = KSPSetType(pcis->ksp_N,KSPPREONLY);CHKERRQ(ierr); 149 ierr = KSPSetFromOptions(pcis->ksp_N);CHKERRQ(ierr); 150 { 151 PetscTruth damp_fixed, 152 remove_nullspace_fixed, 153 set_damping_factor_floating, 154 not_damp_floating, 155 not_remove_nullspace_floating; 156 PetscReal fixed_factor, 157 floating_factor; 158 159 ierr = PetscOptionsGetReal(pc_ctx->prefix,"-pc_is_damp_fixed",&fixed_factor,&damp_fixed);CHKERRQ(ierr); 160 if (!damp_fixed) { fixed_factor = 0.0; } 161 ierr = PetscOptionsHasName(pc_ctx->prefix,"-pc_is_damp_fixed",&damp_fixed);CHKERRQ(ierr); 162 163 ierr = PetscOptionsHasName(pc_ctx->prefix,"-pc_is_remove_nullspace_fixed",&remove_nullspace_fixed);CHKERRQ(ierr); 164 165 ierr = PetscOptionsGetReal(pc_ctx->prefix,"-pc_is_set_damping_factor_floating", 166 &floating_factor,&set_damping_factor_floating);CHKERRQ(ierr); 167 if (!set_damping_factor_floating) { floating_factor = 0.0; } 168 ierr = PetscOptionsHasName(pc_ctx->prefix,"-pc_is_set_damping_factor_floating",&set_damping_factor_floating);CHKERRQ(ierr); 169 if (!set_damping_factor_floating) { floating_factor = 1.e-12; } 170 171 ierr = PetscOptionsHasName(pc_ctx->prefix,"-pc_is_not_damp_floating",¬_damp_floating);CHKERRQ(ierr); 172 173 ierr = PetscOptionsHasName(pc_ctx->prefix,"-pc_is_not_remove_nullspace_floating",¬_remove_nullspace_floating);CHKERRQ(ierr); 174 175 if (pcis->pure_neumann) { /* floating subdomain */ 176 if (!(not_damp_floating)) { 177 ierr = PCLUSetDamping (pc_ctx,floating_factor);CHKERRQ(ierr); 178 ierr = PCILUSetDamping(pc_ctx,floating_factor);CHKERRQ(ierr); 179 } 180 if (!(not_remove_nullspace_floating)){ 181 MatNullSpace nullsp; 182 ierr = MatNullSpaceCreate(PETSC_COMM_SELF,1,0,PETSC_NULL,&nullsp);CHKERRQ(ierr); 183 ierr = PCNullSpaceAttach(pc_ctx,nullsp);CHKERRQ(ierr); 184 ierr = MatNullSpaceDestroy(nullsp);CHKERRQ(ierr); 185 } 186 } else { /* fixed subdomain */ 187 if (damp_fixed) { 188 ierr = PCLUSetDamping (pc_ctx,fixed_factor);CHKERRQ(ierr); 189 ierr = PCILUSetDamping(pc_ctx,fixed_factor);CHKERRQ(ierr); 190 } 191 if (remove_nullspace_fixed) { 192 MatNullSpace nullsp; 193 ierr = MatNullSpaceCreate(PETSC_COMM_SELF,1,0,PETSC_NULL,&nullsp);CHKERRQ(ierr); 194 ierr = PCNullSpaceAttach(pc_ctx,nullsp);CHKERRQ(ierr); 195 ierr = MatNullSpaceDestroy(nullsp);CHKERRQ(ierr); 196 } 197 } 198 } 199 /* the vectors in the following line are dummy arguments, just telling the KSP the vector size. Values are not used */ 200 ierr = KSPSetRhs(pcis->ksp_N,pcis->vec1_N);CHKERRQ(ierr); 201 ierr = KSPSetSolution(pcis->ksp_N,pcis->vec2_N);CHKERRQ(ierr); 202 ierr = KSPSetUp(pcis->ksp_N);CHKERRQ(ierr); 203 } 204 205 ierr = ISLocalToGlobalMappingGetInfo(((Mat_IS*)(pc->mat->data))->mapping,&(pcis->n_neigh),&(pcis->neigh), 206 &(pcis->n_shared),&(pcis->shared));CHKERRQ(ierr); 207 pcis->ISLocalToGlobalMappingGetInfoWasCalled = PETSC_TRUE; 208 209 PetscFunctionReturn(0); 210 } 211 212 /* -------------------------------------------------------------------------- */ 213 /* 214 PCISDestroy - 215 */ 216 #undef __FUNCT__ 217 #define __FUNCT__ "PCISDestroy" 218 int PCISDestroy(PC pc) 219 { 220 PC_IS *pcis = (PC_IS*)(pc->data); 221 int ierr; 222 223 PetscFunctionBegin; 224 225 if (pcis->is_B_local) {ierr = ISDestroy(pcis->is_B_local);CHKERRQ(ierr);} 226 if (pcis->is_I_local) {ierr = ISDestroy(pcis->is_I_local);CHKERRQ(ierr);} 227 if (pcis->is_B_global) {ierr = ISDestroy(pcis->is_B_global);CHKERRQ(ierr);} 228 if (pcis->is_I_global) {ierr = ISDestroy(pcis->is_I_global);CHKERRQ(ierr);} 229 if (pcis->A_II) {ierr = MatDestroy(pcis->A_II);CHKERRQ(ierr);} 230 if (pcis->A_IB) {ierr = MatDestroy(pcis->A_IB);CHKERRQ(ierr);} 231 if (pcis->A_BI) {ierr = MatDestroy(pcis->A_BI);CHKERRQ(ierr);} 232 if (pcis->A_BB) {ierr = MatDestroy(pcis->A_BB);CHKERRQ(ierr);} 233 if (pcis->D) {ierr = VecDestroy(pcis->D);CHKERRQ(ierr);} 234 if (pcis->ksp_N) {ierr = KSPDestroy(pcis->ksp_N);CHKERRQ(ierr);} 235 if (pcis->ksp_D) {ierr = KSPDestroy(pcis->ksp_D);CHKERRQ(ierr);} 236 if (pcis->vec1_N) {ierr = VecDestroy(pcis->vec1_N);CHKERRQ(ierr);} 237 if (pcis->vec2_N) {ierr = VecDestroy(pcis->vec2_N);CHKERRQ(ierr);} 238 if (pcis->vec1_D) {ierr = VecDestroy(pcis->vec1_D);CHKERRQ(ierr);} 239 if (pcis->vec2_D) {ierr = VecDestroy(pcis->vec2_D);CHKERRQ(ierr);} 240 if (pcis->vec3_D) {ierr = VecDestroy(pcis->vec3_D);CHKERRQ(ierr);} 241 if (pcis->vec1_B) {ierr = VecDestroy(pcis->vec1_B);CHKERRQ(ierr);} 242 if (pcis->vec2_B) {ierr = VecDestroy(pcis->vec2_B);CHKERRQ(ierr);} 243 if (pcis->vec3_B) {ierr = VecDestroy(pcis->vec3_B);CHKERRQ(ierr);} 244 if (pcis->vec1_global) {ierr = VecDestroy(pcis->vec1_global);CHKERRQ(ierr);} 245 if (pcis->work_N) {ierr = PetscFree(pcis->work_N);CHKERRQ(ierr);} 246 if (pcis->global_to_D) {ierr = VecScatterDestroy(pcis->global_to_D);CHKERRQ(ierr);} 247 if (pcis->N_to_B) {ierr = VecScatterDestroy(pcis->N_to_B);CHKERRQ(ierr);} 248 if (pcis->global_to_B) {ierr = VecScatterDestroy(pcis->global_to_B);CHKERRQ(ierr);} 249 if (pcis->ISLocalToGlobalMappingGetInfoWasCalled) { 250 ierr = ISLocalToGlobalMappingRestoreInfo((ISLocalToGlobalMapping)0,&(pcis->n_neigh),&(pcis->neigh),&(pcis->n_shared),&(pcis->shared));CHKERRQ(ierr); 251 } 252 253 PetscFunctionReturn(0); 254 } 255 256 /* -------------------------------------------------------------------------- */ 257 /* 258 PCISCreate - 259 */ 260 #undef __FUNCT__ 261 #define __FUNCT__ "PCISCreate" 262 int PCISCreate(PC pc) 263 { 264 PC_IS *pcis = (PC_IS*)(pc->data); 265 266 PetscFunctionBegin; 267 268 pcis->is_B_local = 0; 269 pcis->is_I_local = 0; 270 pcis->is_B_global = 0; 271 pcis->is_I_global = 0; 272 pcis->A_II = 0; 273 pcis->A_IB = 0; 274 pcis->A_BI = 0; 275 pcis->A_BB = 0; 276 pcis->D = 0; 277 pcis->ksp_N = 0; 278 pcis->ksp_D = 0; 279 pcis->vec1_N = 0; 280 pcis->vec2_N = 0; 281 pcis->vec1_D = 0; 282 pcis->vec2_D = 0; 283 pcis->vec3_D = 0; 284 pcis->vec1_B = 0; 285 pcis->vec2_B = 0; 286 pcis->vec3_B = 0; 287 pcis->vec1_global = 0; 288 pcis->work_N = 0; 289 pcis->global_to_D = 0; 290 pcis->N_to_B = 0; 291 pcis->global_to_B = 0; 292 pcis->ISLocalToGlobalMappingGetInfoWasCalled = PETSC_FALSE; 293 294 PetscFunctionReturn(0); 295 } 296 297 /* -------------------------------------------------------------------------- */ 298 /* 299 PCISApplySchur - 300 301 Input parameters: 302 . pc - preconditioner context 303 . v - vector to which the Schur complement is to be applied (it is NOT modified inside this function, UNLESS vec2_B is null) 304 305 Output parameters: 306 . vec1_B - result of Schur complement applied to chunk 307 . vec2_B - garbage (used as work space), or null (and v is used as workspace) 308 . vec1_D - garbage (used as work space) 309 . vec2_D - garbage (used as work space) 310 311 */ 312 #undef __FUNCT__ 313 #define __FUNCT__ "PCIterSuApplySchur" 314 int PCISApplySchur(PC pc, Vec v, Vec vec1_B, Vec vec2_B, Vec vec1_D, Vec vec2_D) 315 { 316 int ierr; 317 PetscScalar m_one = -1.0; 318 PC_IS *pcis = (PC_IS*)(pc->data); 319 320 PetscFunctionBegin; 321 322 if (vec2_B == (Vec)0) { vec2_B = v; } 323 324 ierr = MatMult(pcis->A_BB,v,vec1_B);CHKERRQ(ierr); 325 ierr = MatMult(pcis->A_IB,v,vec1_D);CHKERRQ(ierr); 326 ierr = KSPSetRhs(pcis->ksp_D,vec1_D);CHKERRQ(ierr); 327 ierr = KSPSetSolution(pcis->ksp_D,vec2_D);CHKERRQ(ierr); 328 ierr = KSPSolve(pcis->ksp_D);CHKERRQ(ierr); 329 ierr = MatMult(pcis->A_BI,vec2_D,vec2_B);CHKERRQ(ierr); 330 ierr = VecAXPY(&m_one,vec2_B,vec1_B);CHKERRQ(ierr); 331 332 PetscFunctionReturn(0); 333 } 334 335 /* -------------------------------------------------------------------------- */ 336 /* 337 PCISScatterArrayNToVecB - Scatters interface node values from a big array (of all local nodes, interior or interface, 338 including ghosts) into an interface vector, when in SCATTER_FORWARD mode, or vice-versa, when in SCATTER_REVERSE 339 mode. 340 341 Input parameters: 342 . pc - preconditioner context 343 . array_N - [when in SCATTER_FORWARD mode] Array to be scattered into the vector 344 . v_B - [when in SCATTER_REVERSE mode] Vector to be scattered into the array 345 346 Output parameter: 347 . array_N - [when in SCATTER_REVERSE mode] Array to receive the scattered vector 348 . v_B - [when in SCATTER_FORWARD mode] Vector to receive the scattered array 349 350 Notes: 351 The entries in the array that do not correspond to interface nodes remain unaltered. 352 */ 353 #undef __FUNCT__ 354 #define __FUNCT__ "PCISScatterArrayNToVecB" 355 int PCISScatterArrayNToVecB (PetscScalar *array_N, Vec v_B, InsertMode imode, ScatterMode smode, PC pc) 356 { 357 int i, ierr, *idex; 358 PetscScalar *array_B; 359 PC_IS *pcis = (PC_IS*)(pc->data); 360 361 PetscFunctionBegin; 362 363 ierr = VecGetArray(v_B,&array_B);CHKERRQ(ierr); 364 ierr = ISGetIndices(pcis->is_B_local,&idex);CHKERRQ(ierr); 365 366 if (smode == SCATTER_FORWARD) { 367 if (imode == INSERT_VALUES) { 368 for (i=0; i<pcis->n_B; i++) { array_B[i] = array_N[idex[i]]; } 369 } else { /* ADD_VALUES */ 370 for (i=0; i<pcis->n_B; i++) { array_B[i] += array_N[idex[i]]; } 371 } 372 } else { /* SCATTER_REVERSE */ 373 if (imode == INSERT_VALUES) { 374 for (i=0; i<pcis->n_B; i++) { array_N[idex[i]] = array_B[i]; } 375 } else { /* ADD_VALUES */ 376 for (i=0; i<pcis->n_B; i++) { array_N[idex[i]] += array_B[i]; } 377 } 378 } 379 380 ierr = ISRestoreIndices(pcis->is_B_local,&idex);CHKERRQ(ierr); 381 ierr = VecRestoreArray(v_B,&array_B);CHKERRQ(ierr); 382 383 PetscFunctionReturn(0); 384 } 385 386 /* -------------------------------------------------------------------------- */ 387 /* 388 PCISApplyInvSchur - Solves the Neumann problem related to applying the inverse of the Schur complement. 389 More precisely, solves the problem: 390 [ A_II A_IB ] [ . ] [ 0 ] 391 [ ] [ ] = [ ] 392 [ A_BI A_BB ] [ x ] [ b ] 393 394 Input parameters: 395 . pc - preconditioner context 396 . b - vector of local interface nodes (including ghosts) 397 398 Output parameters: 399 . x - vector of local interface nodes (including ghosts); returns the application of the inverse of the Schur 400 complement to b 401 . vec1_N - vector of local nodes (interior and interface, including ghosts); returns garbage (used as work space) 402 . vec2_N - vector of local nodes (interior and interface, including ghosts); returns garbage (used as work space) 403 404 */ 405 #undef __FUNCT__ 406 #define __FUNCT__ "PCISApplyInvSchur" 407 int PCISApplyInvSchur (PC pc, Vec b, Vec x, Vec vec1_N, Vec vec2_N) 408 { 409 int ierr; 410 PC_IS *pcis = (PC_IS*)(pc->data); 411 PetscScalar zero = 0.0; 412 413 PetscFunctionBegin; 414 415 /* 416 Neumann solvers. 417 Applying the inverse of the local Schur complement, i.e, solving a Neumann 418 Problem with zero at the interior nodes of the RHS and extracting the interface 419 part of the solution. inverse Schur complement is applied to b and the result 420 is stored in x. 421 */ 422 /* Setting the RHS vec1_N */ 423 ierr = VecSet(&zero,vec1_N);CHKERRQ(ierr); 424 ierr = VecScatterBegin(b,vec1_N,INSERT_VALUES,SCATTER_REVERSE,pcis->N_to_B);CHKERRQ(ierr); 425 ierr = VecScatterEnd (b,vec1_N,INSERT_VALUES,SCATTER_REVERSE,pcis->N_to_B);CHKERRQ(ierr); 426 /* Checking for consistency of the RHS */ 427 { 428 PetscTruth flg; 429 ierr = PetscOptionsHasName(PETSC_NULL,"-pc_is_check_consistency",&flg);CHKERRQ(ierr); 430 if (flg) { 431 PetscScalar average; 432 ierr = VecSum(vec1_N,&average);CHKERRQ(ierr); 433 average = average / ((PetscReal)pcis->n); 434 if (pcis->pure_neumann) { 435 ierr = PetscViewerASCIISynchronizedPrintf(PETSC_VIEWER_STDOUT_(pc->comm),"Subdomain %04d is floating. Average = % 1.14e\n", 436 PetscGlobalRank,PetscAbsScalar(average));CHKERRQ(ierr); 437 } else { 438 ierr = PetscViewerASCIISynchronizedPrintf(PETSC_VIEWER_STDOUT_(pc->comm),"Subdomain %04d is fixed. Average = % 1.14e\n", 439 PetscGlobalRank,PetscAbsScalar(average));CHKERRQ(ierr); 440 } 441 PetscViewerFlush(PETSC_VIEWER_STDOUT_(pc->comm)); 442 } 443 } 444 /* Solving the system for vec2_N */ 445 ierr = KSPSetRhs(pcis->ksp_N,vec1_N);CHKERRQ(ierr); 446 ierr = KSPSetSolution(pcis->ksp_N,vec2_N);CHKERRQ(ierr); 447 ierr = KSPSolve(pcis->ksp_N);CHKERRQ(ierr); 448 /* Extracting the local interface vector out of the solution */ 449 ierr = VecScatterBegin(vec2_N,x,INSERT_VALUES,SCATTER_FORWARD,pcis->N_to_B);CHKERRQ(ierr); 450 ierr = VecScatterEnd (vec2_N,x,INSERT_VALUES,SCATTER_FORWARD,pcis->N_to_B);CHKERRQ(ierr); 451 452 PetscFunctionReturn(0); 453 } 454 455 456 457 458 459 460 461 462 463