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