1 /* TODOLIST 2 3 ConstraintsSetup 4 - tolerances for constraints as an option (take care of single precision!) 5 6 Solvers 7 - Add support for cholesky for coarse solver (similar to local solvers) 8 - Propagate ksp prefixes for solvers to mat objects? 9 - Propagate nearnullspace info among levels 10 11 User interface 12 - ** DofSplitting and DM attached to pc? 13 14 Debugging output 15 - * Better management of verbosity levels of debugging output 16 17 Build 18 - make runexe59 19 20 Extra 21 - ** GetRid of PCBDDCApplySchur, use MatSchur instead 22 - *** Is it possible to work with PCBDDCGraph on boundary indices only (less memory consumed)? 23 - add support for computing h,H and related using coordinates? 24 - Change of basis approach does not work with my nonlinear mechanics example. why? (seems not an issue with l2gmap) 25 - Better management in PCIS code 26 - BDDC with MG framework? 27 28 FETIDP 29 - Move FETIDP code to its own classes 30 31 MATIS related operations contained in BDDC code 32 - Provide general case for subassembling 33 - *** Preallocation routines in MatISGetMPIAXAIJ 34 35 */ 36 37 /* ---------------------------------------------------------------------------------------------------------------------------------------------- 38 Implementation of BDDC preconditioner based on: 39 C. Dohrmann "An approximate BDDC preconditioner", Numerical Linear Algebra with Applications Volume 14, Issue 2, pages 149-168, March 2007 40 ---------------------------------------------------------------------------------------------------------------------------------------------- */ 41 42 #include <../src/ksp/pc/impls/bddc/bddc.h> /*I "petscpc.h" I*/ /* includes for fortran wrappers */ 43 #include <../src/ksp/pc/impls/bddc/bddcprivate.h> 44 #include <petscblaslapack.h> 45 46 /* -------------------------------------------------------------------------- */ 47 #undef __FUNCT__ 48 #define __FUNCT__ "PCSetFromOptions_BDDC" 49 PetscErrorCode PCSetFromOptions_BDDC(PC pc) 50 { 51 PC_BDDC *pcbddc = (PC_BDDC*)pc->data; 52 PetscErrorCode ierr; 53 54 PetscFunctionBegin; 55 ierr = PetscOptionsHead("BDDC options");CHKERRQ(ierr); 56 /* Verbose debugging */ 57 ierr = PetscOptionsInt("-pc_bddc_check_level","Verbose output for PCBDDC (intended for debug)","none",pcbddc->dbg_flag,&pcbddc->dbg_flag,NULL);CHKERRQ(ierr); 58 /* Primal space cumstomization */ 59 ierr = PetscOptionsBool("-pc_bddc_use_local_mat_graph","Use or not adjacency graph of local mat for interface analysis","none",pcbddc->use_local_adj,&pcbddc->use_local_adj,NULL);CHKERRQ(ierr); 60 ierr = PetscOptionsBool("-pc_bddc_use_vertices","Use or not corner dofs in coarse space","none",pcbddc->use_vertices,&pcbddc->use_vertices,NULL);CHKERRQ(ierr); 61 ierr = PetscOptionsBool("-pc_bddc_use_edges","Use or not edge constraints in coarse space","none",pcbddc->use_edges,&pcbddc->use_edges,NULL);CHKERRQ(ierr); 62 ierr = PetscOptionsBool("-pc_bddc_use_faces","Use or not face constraints in coarse space","none",pcbddc->use_faces,&pcbddc->use_faces,NULL);CHKERRQ(ierr); 63 /* Change of basis */ 64 ierr = PetscOptionsBool("-pc_bddc_use_change_of_basis","Use or not internal change of basis on local edge nodes","none",pcbddc->use_change_of_basis,&pcbddc->use_change_of_basis,NULL);CHKERRQ(ierr); 65 ierr = PetscOptionsBool("-pc_bddc_use_change_on_faces","Use or not internal change of basis on local face nodes","none",pcbddc->use_change_on_faces,&pcbddc->use_change_on_faces,NULL);CHKERRQ(ierr); 66 if (!pcbddc->use_change_of_basis) { 67 pcbddc->use_change_on_faces = PETSC_FALSE; 68 } 69 /* Switch between M_2 (default) and M_3 preconditioners (as defined by C. Dohrmann in the ref. article) */ 70 ierr = PetscOptionsBool("-pc_bddc_switch_static","Switch on static condensation ops around the interface preconditioner","none",pcbddc->switch_static,&pcbddc->switch_static,NULL);CHKERRQ(ierr); 71 ierr = PetscOptionsInt("-pc_bddc_coarsening_ratio","Set coarsening ratio used in multilevel coarsening","none",pcbddc->coarsening_ratio,&pcbddc->coarsening_ratio,NULL);CHKERRQ(ierr); 72 ierr = PetscOptionsInt("-pc_bddc_levels","Set maximum number of levels for multilevel","none",pcbddc->max_levels,&pcbddc->max_levels,NULL);CHKERRQ(ierr); 73 ierr = PetscOptionsBool("-pc_bddc_use_coarse_estimates","Use estimated eigenvalues for coarse problem","none",pcbddc->use_coarse_estimates,&pcbddc->use_coarse_estimates,NULL);CHKERRQ(ierr); 74 ierr = PetscOptionsBool("-pc_bddc_use_deluxe_scaling","Use deluxe scaling for BDDC","none",pcbddc->use_deluxe_scaling,&pcbddc->use_deluxe_scaling,NULL);CHKERRQ(ierr); 75 ierr = PetscOptionsInt("-pc_bddc_deluxe_threshold","Deluxe subproblems (Schur principal minors) smaller than this value are explicilty computed (-1 computes all)","none",pcbddc->deluxe_threshold,&pcbddc->deluxe_threshold,NULL);CHKERRQ(ierr); 76 ierr = PetscOptionsBool("-pc_bddc_deluxe_rebuild","Whether or not the interface graph for deluxe has to be rebuilt (i.e. use a standard definition of the interface)","none",pcbddc->deluxe_rebuild,&pcbddc->deluxe_rebuild,NULL);CHKERRQ(ierr); 77 ierr = PetscOptionsInt("-pc_bddc_deluxe_layers","Number of dofs' layers for the application of deluxe cheap version (i.e. -1 uses all dofs)","none",pcbddc->deluxe_layers,&pcbddc->deluxe_layers,NULL);CHKERRQ(ierr); 78 ierr = PetscOptionsBool("-pc_bddc_deluxe_use_useradj","Whether or not the CSR graph specified by the user should be used for computing layers (default is to use adj of local mat)","none",pcbddc->deluxe_use_useradj,&pcbddc->deluxe_use_useradj,NULL);CHKERRQ(ierr); 79 ierr = PetscOptionsTail();CHKERRQ(ierr); 80 PetscFunctionReturn(0); 81 } 82 /* -------------------------------------------------------------------------- */ 83 #undef __FUNCT__ 84 #define __FUNCT__ "PCBDDCSetChangeOfBasisLocalMat_BDDC" 85 static PetscErrorCode PCBDDCSetChangeOfBasisLocalMat_BDDC(PC pc, Mat change) 86 { 87 PC_BDDC *pcbddc = (PC_BDDC*)pc->data; 88 PetscErrorCode ierr; 89 90 PetscFunctionBegin; 91 ierr = MatDestroy(&pcbddc->user_ChangeOfBasisMatrix);CHKERRQ(ierr); 92 ierr = PetscObjectReference((PetscObject)change);CHKERRQ(ierr); 93 pcbddc->user_ChangeOfBasisMatrix = change; 94 PetscFunctionReturn(0); 95 } 96 #undef __FUNCT__ 97 #define __FUNCT__ "PCBDDCSetChangeOfBasisLocalMat" 98 /*@ 99 PCBDDCSetChangeOfBasisLocalMat - Set user defined change of basis for local boundary dofs 100 101 Collective on PC 102 103 Input Parameters: 104 + pc - the preconditioning context 105 - change - the local change of basis matrix, either in local (internal + boundary) or in local boundary numbering 106 107 Level: intermediate 108 109 Notes: 110 111 .seealso: PCBDDC 112 @*/ 113 PetscErrorCode PCBDDCSetChangeOfBasisLocalMat(PC pc, Mat change) 114 { 115 PetscErrorCode ierr; 116 117 PetscFunctionBegin; 118 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 119 PetscValidHeaderSpecific(change,MAT_CLASSID,2); 120 ierr = PetscTryMethod(pc,"PCBDDCSetChangeOfBasisLocalMat_C",(PC,Mat),(pc,change));CHKERRQ(ierr); 121 PetscFunctionReturn(0); 122 } 123 /* -------------------------------------------------------------------------- */ 124 #undef __FUNCT__ 125 #define __FUNCT__ "PCBDDCSetPrimalVerticesLocalIS_BDDC" 126 static PetscErrorCode PCBDDCSetPrimalVerticesLocalIS_BDDC(PC pc, IS PrimalVertices) 127 { 128 PC_BDDC *pcbddc = (PC_BDDC*)pc->data; 129 PetscErrorCode ierr; 130 131 PetscFunctionBegin; 132 ierr = ISDestroy(&pcbddc->user_primal_vertices);CHKERRQ(ierr); 133 ierr = PetscObjectReference((PetscObject)PrimalVertices);CHKERRQ(ierr); 134 pcbddc->user_primal_vertices = PrimalVertices; 135 PetscFunctionReturn(0); 136 } 137 #undef __FUNCT__ 138 #define __FUNCT__ "PCBDDCSetPrimalVerticesLocalIS" 139 /*@ 140 PCBDDCSetPrimalVerticesLocalIS - Set additional user defined primal vertices in PCBDDC 141 142 Not collective 143 144 Input Parameters: 145 + pc - the preconditioning context 146 - PrimalVertices - index set of primal vertices in local numbering 147 148 Level: intermediate 149 150 Notes: 151 152 .seealso: PCBDDC 153 @*/ 154 PetscErrorCode PCBDDCSetPrimalVerticesLocalIS(PC pc, IS PrimalVertices) 155 { 156 PetscErrorCode ierr; 157 158 PetscFunctionBegin; 159 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 160 PetscValidHeaderSpecific(PrimalVertices,IS_CLASSID,2); 161 ierr = PetscTryMethod(pc,"PCBDDCSetPrimalVerticesLocalIS_C",(PC,IS),(pc,PrimalVertices));CHKERRQ(ierr); 162 PetscFunctionReturn(0); 163 } 164 /* -------------------------------------------------------------------------- */ 165 #undef __FUNCT__ 166 #define __FUNCT__ "PCBDDCSetCoarseningRatio_BDDC" 167 static PetscErrorCode PCBDDCSetCoarseningRatio_BDDC(PC pc,PetscInt k) 168 { 169 PC_BDDC *pcbddc = (PC_BDDC*)pc->data; 170 171 PetscFunctionBegin; 172 pcbddc->coarsening_ratio = k; 173 PetscFunctionReturn(0); 174 } 175 176 #undef __FUNCT__ 177 #define __FUNCT__ "PCBDDCSetCoarseningRatio" 178 /*@ 179 PCBDDCSetCoarseningRatio - Set coarsening ratio used in multilevel 180 181 Logically collective on PC 182 183 Input Parameters: 184 + pc - the preconditioning context 185 - k - coarsening ratio (H/h at the coarser level) 186 187 Approximatively k subdomains at the finer level will be aggregated into a single subdomain at the coarser level 188 189 Level: intermediate 190 191 Notes: 192 193 .seealso: PCBDDC 194 @*/ 195 PetscErrorCode PCBDDCSetCoarseningRatio(PC pc,PetscInt k) 196 { 197 PetscErrorCode ierr; 198 199 PetscFunctionBegin; 200 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 201 PetscValidLogicalCollectiveInt(pc,k,2); 202 ierr = PetscTryMethod(pc,"PCBDDCSetCoarseningRatio_C",(PC,PetscInt),(pc,k));CHKERRQ(ierr); 203 PetscFunctionReturn(0); 204 } 205 206 /* The following functions (PCBDDCSetUseExactDirichlet PCBDDCSetLevel) are not public */ 207 #undef __FUNCT__ 208 #define __FUNCT__ "PCBDDCSetUseExactDirichlet_BDDC" 209 static PetscErrorCode PCBDDCSetUseExactDirichlet_BDDC(PC pc,PetscBool flg) 210 { 211 PC_BDDC *pcbddc = (PC_BDDC*)pc->data; 212 213 PetscFunctionBegin; 214 pcbddc->use_exact_dirichlet_trick = flg; 215 PetscFunctionReturn(0); 216 } 217 218 #undef __FUNCT__ 219 #define __FUNCT__ "PCBDDCSetUseExactDirichlet" 220 PetscErrorCode PCBDDCSetUseExactDirichlet(PC pc,PetscBool flg) 221 { 222 PetscErrorCode ierr; 223 224 PetscFunctionBegin; 225 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 226 PetscValidLogicalCollectiveBool(pc,flg,2); 227 ierr = PetscTryMethod(pc,"PCBDDCSetUseExactDirichlet_C",(PC,PetscBool),(pc,flg));CHKERRQ(ierr); 228 PetscFunctionReturn(0); 229 } 230 231 #undef __FUNCT__ 232 #define __FUNCT__ "PCBDDCSetLevel_BDDC" 233 static PetscErrorCode PCBDDCSetLevel_BDDC(PC pc,PetscInt level) 234 { 235 PC_BDDC *pcbddc = (PC_BDDC*)pc->data; 236 237 PetscFunctionBegin; 238 pcbddc->current_level = level; 239 PetscFunctionReturn(0); 240 } 241 242 #undef __FUNCT__ 243 #define __FUNCT__ "PCBDDCSetLevel" 244 PetscErrorCode PCBDDCSetLevel(PC pc,PetscInt level) 245 { 246 PetscErrorCode ierr; 247 248 PetscFunctionBegin; 249 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 250 PetscValidLogicalCollectiveInt(pc,level,2); 251 ierr = PetscTryMethod(pc,"PCBDDCSetLevel_C",(PC,PetscInt),(pc,level));CHKERRQ(ierr); 252 PetscFunctionReturn(0); 253 } 254 255 #undef __FUNCT__ 256 #define __FUNCT__ "PCBDDCSetLevels_BDDC" 257 static PetscErrorCode PCBDDCSetLevels_BDDC(PC pc,PetscInt levels) 258 { 259 PC_BDDC *pcbddc = (PC_BDDC*)pc->data; 260 261 PetscFunctionBegin; 262 pcbddc->max_levels = levels; 263 PetscFunctionReturn(0); 264 } 265 266 #undef __FUNCT__ 267 #define __FUNCT__ "PCBDDCSetLevels" 268 /*@ 269 PCBDDCSetLevels - Sets the maximum number of levels for multilevel 270 271 Logically collective on PC 272 273 Input Parameters: 274 + pc - the preconditioning context 275 - levels - the maximum number of levels (max 9) 276 277 Default value is 0, i.e. traditional one-level BDDC 278 279 Level: intermediate 280 281 Notes: 282 283 .seealso: PCBDDC 284 @*/ 285 PetscErrorCode PCBDDCSetLevels(PC pc,PetscInt levels) 286 { 287 PetscErrorCode ierr; 288 289 PetscFunctionBegin; 290 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 291 PetscValidLogicalCollectiveInt(pc,levels,2); 292 if (levels > 99) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"Maximum number of levels for bddc is 99\n"); 293 ierr = PetscTryMethod(pc,"PCBDDCSetLevels_C",(PC,PetscInt),(pc,levels));CHKERRQ(ierr); 294 PetscFunctionReturn(0); 295 } 296 /* -------------------------------------------------------------------------- */ 297 298 #undef __FUNCT__ 299 #define __FUNCT__ "PCBDDCSetNullSpace_BDDC" 300 static PetscErrorCode PCBDDCSetNullSpace_BDDC(PC pc,MatNullSpace NullSpace) 301 { 302 PC_BDDC *pcbddc = (PC_BDDC*)pc->data; 303 PetscErrorCode ierr; 304 305 PetscFunctionBegin; 306 ierr = PetscObjectReference((PetscObject)NullSpace);CHKERRQ(ierr); 307 ierr = MatNullSpaceDestroy(&pcbddc->NullSpace);CHKERRQ(ierr); 308 pcbddc->NullSpace = NullSpace; 309 PetscFunctionReturn(0); 310 } 311 312 #undef __FUNCT__ 313 #define __FUNCT__ "PCBDDCSetNullSpace" 314 /*@ 315 PCBDDCSetNullSpace - Set nullspace for BDDC operator 316 317 Logically collective on PC and MatNullSpace 318 319 Input Parameters: 320 + pc - the preconditioning context 321 - NullSpace - Null space of the linear operator to be preconditioned (Pmat) 322 323 Level: intermediate 324 325 Notes: 326 327 .seealso: PCBDDC 328 @*/ 329 PetscErrorCode PCBDDCSetNullSpace(PC pc,MatNullSpace NullSpace) 330 { 331 PetscErrorCode ierr; 332 333 PetscFunctionBegin; 334 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 335 PetscValidHeaderSpecific(NullSpace,MAT_NULLSPACE_CLASSID,2); 336 PetscCheckSameComm(pc,1,NullSpace,2); 337 ierr = PetscTryMethod(pc,"PCBDDCSetNullSpace_C",(PC,MatNullSpace),(pc,NullSpace));CHKERRQ(ierr); 338 PetscFunctionReturn(0); 339 } 340 /* -------------------------------------------------------------------------- */ 341 342 #undef __FUNCT__ 343 #define __FUNCT__ "PCBDDCSetDirichletBoundaries_BDDC" 344 static PetscErrorCode PCBDDCSetDirichletBoundaries_BDDC(PC pc,IS DirichletBoundaries) 345 { 346 PC_BDDC *pcbddc = (PC_BDDC*)pc->data; 347 PetscErrorCode ierr; 348 349 PetscFunctionBegin; 350 /* last user setting takes precendence -> destroy any other customization */ 351 ierr = ISDestroy(&pcbddc->DirichletBoundariesLocal);CHKERRQ(ierr); 352 ierr = ISDestroy(&pcbddc->DirichletBoundaries);CHKERRQ(ierr); 353 ierr = PetscObjectReference((PetscObject)DirichletBoundaries);CHKERRQ(ierr); 354 pcbddc->DirichletBoundaries = DirichletBoundaries; 355 pcbddc->recompute_topography = PETSC_TRUE; 356 PetscFunctionReturn(0); 357 } 358 359 #undef __FUNCT__ 360 #define __FUNCT__ "PCBDDCSetDirichletBoundaries" 361 /*@ 362 PCBDDCSetDirichletBoundaries - Set IS defining Dirichlet boundaries for the global problem. 363 364 Collective 365 366 Input Parameters: 367 + pc - the preconditioning context 368 - DirichletBoundaries - parallel IS defining the Dirichlet boundaries 369 370 Level: intermediate 371 372 Notes: Any process can list any global node 373 374 .seealso: PCBDDC 375 @*/ 376 PetscErrorCode PCBDDCSetDirichletBoundaries(PC pc,IS DirichletBoundaries) 377 { 378 PetscErrorCode ierr; 379 380 PetscFunctionBegin; 381 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 382 PetscValidHeaderSpecific(DirichletBoundaries,IS_CLASSID,2); 383 PetscCheckSameComm(pc,1,DirichletBoundaries,2); 384 ierr = PetscTryMethod(pc,"PCBDDCSetDirichletBoundaries_C",(PC,IS),(pc,DirichletBoundaries));CHKERRQ(ierr); 385 PetscFunctionReturn(0); 386 } 387 /* -------------------------------------------------------------------------- */ 388 389 #undef __FUNCT__ 390 #define __FUNCT__ "PCBDDCSetDirichletBoundariesLocal_BDDC" 391 static PetscErrorCode PCBDDCSetDirichletBoundariesLocal_BDDC(PC pc,IS DirichletBoundaries) 392 { 393 PC_BDDC *pcbddc = (PC_BDDC*)pc->data; 394 PetscErrorCode ierr; 395 396 PetscFunctionBegin; 397 /* last user setting takes precendence -> destroy any other customization */ 398 ierr = ISDestroy(&pcbddc->DirichletBoundariesLocal);CHKERRQ(ierr); 399 ierr = ISDestroy(&pcbddc->DirichletBoundaries);CHKERRQ(ierr); 400 ierr = PetscObjectReference((PetscObject)DirichletBoundaries);CHKERRQ(ierr); 401 pcbddc->DirichletBoundariesLocal = DirichletBoundaries; 402 pcbddc->recompute_topography = PETSC_TRUE; 403 PetscFunctionReturn(0); 404 } 405 406 #undef __FUNCT__ 407 #define __FUNCT__ "PCBDDCSetDirichletBoundariesLocal" 408 /*@ 409 PCBDDCSetDirichletBoundariesLocal - Set IS defining Dirichlet boundaries for the global problem in local ordering. 410 411 Collective 412 413 Input Parameters: 414 + pc - the preconditioning context 415 - DirichletBoundaries - parallel IS defining the Dirichlet boundaries (in local ordering) 416 417 Level: intermediate 418 419 Notes: 420 421 .seealso: PCBDDC 422 @*/ 423 PetscErrorCode PCBDDCSetDirichletBoundariesLocal(PC pc,IS DirichletBoundaries) 424 { 425 PetscErrorCode ierr; 426 427 PetscFunctionBegin; 428 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 429 PetscValidHeaderSpecific(DirichletBoundaries,IS_CLASSID,2); 430 PetscCheckSameComm(pc,1,DirichletBoundaries,2); 431 ierr = PetscTryMethod(pc,"PCBDDCSetDirichletBoundariesLocal_C",(PC,IS),(pc,DirichletBoundaries));CHKERRQ(ierr); 432 PetscFunctionReturn(0); 433 } 434 /* -------------------------------------------------------------------------- */ 435 436 #undef __FUNCT__ 437 #define __FUNCT__ "PCBDDCSetNeumannBoundaries_BDDC" 438 static PetscErrorCode PCBDDCSetNeumannBoundaries_BDDC(PC pc,IS NeumannBoundaries) 439 { 440 PC_BDDC *pcbddc = (PC_BDDC*)pc->data; 441 PetscErrorCode ierr; 442 443 PetscFunctionBegin; 444 /* last user setting takes precendence -> destroy any other customization */ 445 ierr = ISDestroy(&pcbddc->NeumannBoundariesLocal);CHKERRQ(ierr); 446 ierr = ISDestroy(&pcbddc->NeumannBoundaries);CHKERRQ(ierr); 447 ierr = PetscObjectReference((PetscObject)NeumannBoundaries);CHKERRQ(ierr); 448 pcbddc->NeumannBoundaries = NeumannBoundaries; 449 pcbddc->recompute_topography = PETSC_TRUE; 450 PetscFunctionReturn(0); 451 } 452 453 #undef __FUNCT__ 454 #define __FUNCT__ "PCBDDCSetNeumannBoundaries" 455 /*@ 456 PCBDDCSetNeumannBoundaries - Set IS defining Neumann boundaries for the global problem. 457 458 Collective 459 460 Input Parameters: 461 + pc - the preconditioning context 462 - NeumannBoundaries - parallel IS defining the Neumann boundaries 463 464 Level: intermediate 465 466 Notes: Any process can list any global node 467 468 .seealso: PCBDDC 469 @*/ 470 PetscErrorCode PCBDDCSetNeumannBoundaries(PC pc,IS NeumannBoundaries) 471 { 472 PetscErrorCode ierr; 473 474 PetscFunctionBegin; 475 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 476 PetscValidHeaderSpecific(NeumannBoundaries,IS_CLASSID,2); 477 PetscCheckSameComm(pc,1,NeumannBoundaries,2); 478 ierr = PetscTryMethod(pc,"PCBDDCSetNeumannBoundaries_C",(PC,IS),(pc,NeumannBoundaries));CHKERRQ(ierr); 479 PetscFunctionReturn(0); 480 } 481 /* -------------------------------------------------------------------------- */ 482 483 #undef __FUNCT__ 484 #define __FUNCT__ "PCBDDCSetNeumannBoundariesLocal_BDDC" 485 static PetscErrorCode PCBDDCSetNeumannBoundariesLocal_BDDC(PC pc,IS NeumannBoundaries) 486 { 487 PC_BDDC *pcbddc = (PC_BDDC*)pc->data; 488 PetscErrorCode ierr; 489 490 PetscFunctionBegin; 491 /* last user setting takes precendence -> destroy any other customization */ 492 ierr = ISDestroy(&pcbddc->NeumannBoundariesLocal);CHKERRQ(ierr); 493 ierr = ISDestroy(&pcbddc->NeumannBoundaries);CHKERRQ(ierr); 494 ierr = PetscObjectReference((PetscObject)NeumannBoundaries);CHKERRQ(ierr); 495 pcbddc->NeumannBoundariesLocal = NeumannBoundaries; 496 pcbddc->recompute_topography = PETSC_TRUE; 497 PetscFunctionReturn(0); 498 } 499 500 #undef __FUNCT__ 501 #define __FUNCT__ "PCBDDCSetNeumannBoundariesLocal" 502 /*@ 503 PCBDDCSetNeumannBoundariesLocal - Set IS defining Neumann boundaries for the global problem in local ordering. 504 505 Collective 506 507 Input Parameters: 508 + pc - the preconditioning context 509 - NeumannBoundaries - parallel IS defining the subdomain part of Neumann boundaries (in local ordering) 510 511 Level: intermediate 512 513 Notes: 514 515 .seealso: PCBDDC 516 @*/ 517 PetscErrorCode PCBDDCSetNeumannBoundariesLocal(PC pc,IS NeumannBoundaries) 518 { 519 PetscErrorCode ierr; 520 521 PetscFunctionBegin; 522 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 523 PetscValidHeaderSpecific(NeumannBoundaries,IS_CLASSID,2); 524 PetscCheckSameComm(pc,1,NeumannBoundaries,2); 525 ierr = PetscTryMethod(pc,"PCBDDCSetNeumannBoundariesLocal_C",(PC,IS),(pc,NeumannBoundaries));CHKERRQ(ierr); 526 PetscFunctionReturn(0); 527 } 528 /* -------------------------------------------------------------------------- */ 529 530 #undef __FUNCT__ 531 #define __FUNCT__ "PCBDDCGetDirichletBoundaries_BDDC" 532 static PetscErrorCode PCBDDCGetDirichletBoundaries_BDDC(PC pc,IS *DirichletBoundaries) 533 { 534 PC_BDDC *pcbddc = (PC_BDDC*)pc->data; 535 536 PetscFunctionBegin; 537 *DirichletBoundaries = pcbddc->DirichletBoundaries; 538 PetscFunctionReturn(0); 539 } 540 541 #undef __FUNCT__ 542 #define __FUNCT__ "PCBDDCGetDirichletBoundaries" 543 /*@ 544 PCBDDCGetDirichletBoundaries - Get parallel IS for Dirichlet boundaries 545 546 Collective 547 548 Input Parameters: 549 . pc - the preconditioning context 550 551 Output Parameters: 552 . DirichletBoundaries - index set defining the Dirichlet boundaries 553 554 Level: intermediate 555 556 Notes: The IS returned (if any) is the same passed in earlier by the user with PCBDDCSetDirichletBoundaries 557 558 .seealso: PCBDDC 559 @*/ 560 PetscErrorCode PCBDDCGetDirichletBoundaries(PC pc,IS *DirichletBoundaries) 561 { 562 PetscErrorCode ierr; 563 564 PetscFunctionBegin; 565 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 566 ierr = PetscUseMethod(pc,"PCBDDCGetDirichletBoundaries_C",(PC,IS*),(pc,DirichletBoundaries));CHKERRQ(ierr); 567 PetscFunctionReturn(0); 568 } 569 /* -------------------------------------------------------------------------- */ 570 571 #undef __FUNCT__ 572 #define __FUNCT__ "PCBDDCGetDirichletBoundariesLocal_BDDC" 573 static PetscErrorCode PCBDDCGetDirichletBoundariesLocal_BDDC(PC pc,IS *DirichletBoundaries) 574 { 575 PC_BDDC *pcbddc = (PC_BDDC*)pc->data; 576 577 PetscFunctionBegin; 578 *DirichletBoundaries = pcbddc->DirichletBoundariesLocal; 579 PetscFunctionReturn(0); 580 } 581 582 #undef __FUNCT__ 583 #define __FUNCT__ "PCBDDCGetDirichletBoundariesLocal" 584 /*@ 585 PCBDDCGetDirichletBoundariesLocal - Get parallel IS for Dirichlet boundaries (in local ordering) 586 587 Collective 588 589 Input Parameters: 590 . pc - the preconditioning context 591 592 Output Parameters: 593 . DirichletBoundaries - index set defining the subdomain part of Dirichlet boundaries 594 595 Level: intermediate 596 597 Notes: 598 599 .seealso: PCBDDC 600 @*/ 601 PetscErrorCode PCBDDCGetDirichletBoundariesLocal(PC pc,IS *DirichletBoundaries) 602 { 603 PetscErrorCode ierr; 604 605 PetscFunctionBegin; 606 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 607 ierr = PetscUseMethod(pc,"PCBDDCGetDirichletBoundariesLocal_C",(PC,IS*),(pc,DirichletBoundaries));CHKERRQ(ierr); 608 PetscFunctionReturn(0); 609 } 610 /* -------------------------------------------------------------------------- */ 611 612 #undef __FUNCT__ 613 #define __FUNCT__ "PCBDDCGetNeumannBoundaries_BDDC" 614 static PetscErrorCode PCBDDCGetNeumannBoundaries_BDDC(PC pc,IS *NeumannBoundaries) 615 { 616 PC_BDDC *pcbddc = (PC_BDDC*)pc->data; 617 618 PetscFunctionBegin; 619 *NeumannBoundaries = pcbddc->NeumannBoundaries; 620 PetscFunctionReturn(0); 621 } 622 623 #undef __FUNCT__ 624 #define __FUNCT__ "PCBDDCGetNeumannBoundaries" 625 /*@ 626 PCBDDCGetNeumannBoundaries - Get parallel IS for Neumann boundaries 627 628 Collective 629 630 Input Parameters: 631 . pc - the preconditioning context 632 633 Output Parameters: 634 . NeumannBoundaries - index set defining the Neumann boundaries 635 636 Level: intermediate 637 638 Notes: The IS returned (if any) is the same passed in earlier by the user with PCBDDCSetNeumannBoundaries 639 640 .seealso: PCBDDC 641 @*/ 642 PetscErrorCode PCBDDCGetNeumannBoundaries(PC pc,IS *NeumannBoundaries) 643 { 644 PetscErrorCode ierr; 645 646 PetscFunctionBegin; 647 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 648 ierr = PetscUseMethod(pc,"PCBDDCGetNeumannBoundaries_C",(PC,IS*),(pc,NeumannBoundaries));CHKERRQ(ierr); 649 PetscFunctionReturn(0); 650 } 651 /* -------------------------------------------------------------------------- */ 652 653 #undef __FUNCT__ 654 #define __FUNCT__ "PCBDDCGetNeumannBoundariesLocal_BDDC" 655 static PetscErrorCode PCBDDCGetNeumannBoundariesLocal_BDDC(PC pc,IS *NeumannBoundaries) 656 { 657 PC_BDDC *pcbddc = (PC_BDDC*)pc->data; 658 659 PetscFunctionBegin; 660 *NeumannBoundaries = pcbddc->NeumannBoundariesLocal; 661 PetscFunctionReturn(0); 662 } 663 664 #undef __FUNCT__ 665 #define __FUNCT__ "PCBDDCGetNeumannBoundariesLocal" 666 /*@ 667 PCBDDCGetNeumannBoundariesLocal - Get parallel IS for Neumann boundaries (in local ordering) 668 669 Collective 670 671 Input Parameters: 672 . pc - the preconditioning context 673 674 Output Parameters: 675 . NeumannBoundaries - index set defining the subdomain part of Neumann boundaries 676 677 Level: intermediate 678 679 Notes: 680 681 .seealso: PCBDDC 682 @*/ 683 PetscErrorCode PCBDDCGetNeumannBoundariesLocal(PC pc,IS *NeumannBoundaries) 684 { 685 PetscErrorCode ierr; 686 687 PetscFunctionBegin; 688 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 689 ierr = PetscUseMethod(pc,"PCBDDCGetNeumannBoundariesLocal_C",(PC,IS*),(pc,NeumannBoundaries));CHKERRQ(ierr); 690 PetscFunctionReturn(0); 691 } 692 /* -------------------------------------------------------------------------- */ 693 694 #undef __FUNCT__ 695 #define __FUNCT__ "PCBDDCSetLocalAdjacencyGraph_BDDC" 696 static PetscErrorCode PCBDDCSetLocalAdjacencyGraph_BDDC(PC pc, PetscInt nvtxs,const PetscInt xadj[],const PetscInt adjncy[], PetscCopyMode copymode) 697 { 698 PC_BDDC *pcbddc = (PC_BDDC*)pc->data; 699 PCBDDCGraph mat_graph = pcbddc->mat_graph; 700 PetscErrorCode ierr; 701 702 PetscFunctionBegin; 703 /* free old CSR */ 704 ierr = PCBDDCGraphResetCSR(mat_graph);CHKERRQ(ierr); 705 /* TODO: PCBDDCGraphSetAdjacency */ 706 /* get CSR into graph structure */ 707 if (copymode == PETSC_COPY_VALUES) { 708 ierr = PetscMalloc1((nvtxs+1),&mat_graph->xadj);CHKERRQ(ierr); 709 ierr = PetscMalloc1(xadj[nvtxs],&mat_graph->adjncy);CHKERRQ(ierr); 710 ierr = PetscMemcpy(mat_graph->xadj,xadj,(nvtxs+1)*sizeof(PetscInt));CHKERRQ(ierr); 711 ierr = PetscMemcpy(mat_graph->adjncy,adjncy,xadj[nvtxs]*sizeof(PetscInt));CHKERRQ(ierr); 712 } else if (copymode == PETSC_OWN_POINTER) { 713 mat_graph->xadj = (PetscInt*)xadj; 714 mat_graph->adjncy = (PetscInt*)adjncy; 715 } 716 mat_graph->nvtxs_csr = nvtxs; 717 PetscFunctionReturn(0); 718 } 719 720 #undef __FUNCT__ 721 #define __FUNCT__ "PCBDDCSetLocalAdjacencyGraph" 722 /*@ 723 PCBDDCSetLocalAdjacencyGraph - Set adjacency structure (CSR graph) of the local Neumann matrix 724 725 Not collective 726 727 Input Parameters: 728 + pc - the preconditioning context 729 . nvtxs - number of local vertices of the graph (i.e., the local size of your problem) 730 . xadj, adjncy - the CSR graph 731 - copymode - either PETSC_COPY_VALUES or PETSC_OWN_POINTER. 732 733 Level: intermediate 734 735 Notes: 736 737 .seealso: PCBDDC,PetscCopyMode 738 @*/ 739 PetscErrorCode PCBDDCSetLocalAdjacencyGraph(PC pc,PetscInt nvtxs,const PetscInt xadj[],const PetscInt adjncy[], PetscCopyMode copymode) 740 { 741 void (*f)(void) = 0; 742 PetscErrorCode ierr; 743 744 PetscFunctionBegin; 745 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 746 PetscValidIntPointer(xadj,3); 747 PetscValidIntPointer(xadj,4); 748 if (copymode != PETSC_COPY_VALUES && copymode != PETSC_OWN_POINTER) { 749 SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_SUP,"Unsupported copy mode %d in %s\n",copymode,__FUNCT__); 750 } 751 ierr = PetscTryMethod(pc,"PCBDDCSetLocalAdjacencyGraph_C",(PC,PetscInt,const PetscInt[],const PetscInt[],PetscCopyMode),(pc,nvtxs,xadj,adjncy,copymode));CHKERRQ(ierr); 752 /* free arrays if PCBDDC is not the PC type */ 753 ierr = PetscObjectQueryFunction((PetscObject)pc,"PCBDDCSetLocalAdjacencyGraph_C",&f);CHKERRQ(ierr); 754 if (!f && copymode == PETSC_OWN_POINTER) { 755 ierr = PetscFree(xadj);CHKERRQ(ierr); 756 ierr = PetscFree(adjncy);CHKERRQ(ierr); 757 } 758 PetscFunctionReturn(0); 759 } 760 /* -------------------------------------------------------------------------- */ 761 762 #undef __FUNCT__ 763 #define __FUNCT__ "PCBDDCSetDofsSplittingLocal_BDDC" 764 static PetscErrorCode PCBDDCSetDofsSplittingLocal_BDDC(PC pc,PetscInt n_is, IS ISForDofs[]) 765 { 766 PC_BDDC *pcbddc = (PC_BDDC*)pc->data; 767 PetscInt i; 768 PetscErrorCode ierr; 769 770 PetscFunctionBegin; 771 /* Destroy ISes if they were already set */ 772 for (i=0;i<pcbddc->n_ISForDofsLocal;i++) { 773 ierr = ISDestroy(&pcbddc->ISForDofsLocal[i]);CHKERRQ(ierr); 774 } 775 ierr = PetscFree(pcbddc->ISForDofsLocal);CHKERRQ(ierr); 776 /* last user setting takes precendence -> destroy any other customization */ 777 for (i=0;i<pcbddc->n_ISForDofs;i++) { 778 ierr = ISDestroy(&pcbddc->ISForDofs[i]);CHKERRQ(ierr); 779 } 780 ierr = PetscFree(pcbddc->ISForDofs);CHKERRQ(ierr); 781 pcbddc->n_ISForDofs = 0; 782 /* allocate space then set */ 783 ierr = PetscMalloc(n_is*sizeof(IS),&pcbddc->ISForDofsLocal);CHKERRQ(ierr); 784 for (i=0;i<n_is;i++) { 785 ierr = PetscObjectReference((PetscObject)ISForDofs[i]);CHKERRQ(ierr); 786 pcbddc->ISForDofsLocal[i]=ISForDofs[i]; 787 } 788 pcbddc->n_ISForDofsLocal=n_is; 789 if (n_is) pcbddc->user_provided_isfordofs = PETSC_TRUE; 790 pcbddc->recompute_topography = PETSC_TRUE; 791 PetscFunctionReturn(0); 792 } 793 794 #undef __FUNCT__ 795 #define __FUNCT__ "PCBDDCSetDofsSplittingLocal" 796 /*@ 797 PCBDDCSetDofsSplittingLocal - Set index sets defining fields of the local subdomain matrix 798 799 Collective 800 801 Input Parameters: 802 + pc - the preconditioning context 803 - n_is - number of index sets defining the fields 804 . ISForDofs - array of IS describing the fields in local ordering 805 806 Level: intermediate 807 808 Notes: n_is should be the same among processes. Not all nodes need to be listed: unlisted nodes will belong to a different field. 809 810 .seealso: PCBDDC 811 @*/ 812 PetscErrorCode PCBDDCSetDofsSplittingLocal(PC pc,PetscInt n_is, IS ISForDofs[]) 813 { 814 PetscInt i; 815 PetscErrorCode ierr; 816 817 PetscFunctionBegin; 818 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 819 PetscValidLogicalCollectiveInt(pc,n_is,2); 820 for (i=0;i<n_is;i++) { 821 PetscCheckSameComm(pc,1,ISForDofs[i],3); 822 PetscValidHeaderSpecific(ISForDofs[i],IS_CLASSID,3); 823 } 824 ierr = PetscTryMethod(pc,"PCBDDCSetDofsSplitting_C",(PC,PetscInt,IS[]),(pc,n_is,ISForDofs));CHKERRQ(ierr); 825 PetscFunctionReturn(0); 826 } 827 /* -------------------------------------------------------------------------- */ 828 829 #undef __FUNCT__ 830 #define __FUNCT__ "PCBDDCSetDofsSplitting_BDDC" 831 static PetscErrorCode PCBDDCSetDofsSplitting_BDDC(PC pc,PetscInt n_is, IS ISForDofs[]) 832 { 833 PC_BDDC *pcbddc = (PC_BDDC*)pc->data; 834 PetscInt i; 835 PetscErrorCode ierr; 836 837 PetscFunctionBegin; 838 /* Destroy ISes if they were already set */ 839 for (i=0;i<pcbddc->n_ISForDofs;i++) { 840 ierr = ISDestroy(&pcbddc->ISForDofs[i]);CHKERRQ(ierr); 841 } 842 ierr = PetscFree(pcbddc->ISForDofs);CHKERRQ(ierr); 843 /* last user setting takes precendence -> destroy any other customization */ 844 for (i=0;i<pcbddc->n_ISForDofsLocal;i++) { 845 ierr = ISDestroy(&pcbddc->ISForDofsLocal[i]);CHKERRQ(ierr); 846 } 847 ierr = PetscFree(pcbddc->ISForDofsLocal);CHKERRQ(ierr); 848 pcbddc->n_ISForDofsLocal = 0; 849 /* allocate space then set */ 850 ierr = PetscMalloc1(n_is,&pcbddc->ISForDofs);CHKERRQ(ierr); 851 for (i=0;i<n_is;i++) { 852 ierr = PetscObjectReference((PetscObject)ISForDofs[i]);CHKERRQ(ierr); 853 pcbddc->ISForDofs[i]=ISForDofs[i]; 854 } 855 pcbddc->n_ISForDofs=n_is; 856 if (n_is) pcbddc->user_provided_isfordofs = PETSC_TRUE; 857 pcbddc->recompute_topography = PETSC_TRUE; 858 PetscFunctionReturn(0); 859 } 860 861 #undef __FUNCT__ 862 #define __FUNCT__ "PCBDDCSetDofsSplitting" 863 /*@ 864 PCBDDCSetDofsSplitting - Set index sets defining fields of the global matrix 865 866 Collective 867 868 Input Parameters: 869 + pc - the preconditioning context 870 - n_is - number of index sets defining the fields 871 . ISForDofs - array of IS describing the fields in global ordering 872 873 Level: intermediate 874 875 Notes: Any process can list any global node. Not all nodes need to be listed: unlisted nodes will belong to a different field. 876 877 .seealso: PCBDDC 878 @*/ 879 PetscErrorCode PCBDDCSetDofsSplitting(PC pc,PetscInt n_is, IS ISForDofs[]) 880 { 881 PetscInt i; 882 PetscErrorCode ierr; 883 884 PetscFunctionBegin; 885 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 886 PetscValidLogicalCollectiveInt(pc,n_is,2); 887 for (i=0;i<n_is;i++) { 888 PetscCheckSameComm(pc,1,ISForDofs[i],3); 889 PetscValidHeaderSpecific(ISForDofs[i],IS_CLASSID,3); 890 } 891 ierr = PetscTryMethod(pc,"PCBDDCSetDofsSplitting_C",(PC,PetscInt,IS[]),(pc,n_is,ISForDofs));CHKERRQ(ierr); 892 PetscFunctionReturn(0); 893 } 894 /* -------------------------------------------------------------------------- */ 895 #undef __FUNCT__ 896 #define __FUNCT__ "PCPreSolve_BDDC" 897 /* -------------------------------------------------------------------------- */ 898 /* 899 PCPreSolve_BDDC - Changes the right hand side and (if necessary) the initial 900 guess if a transformation of basis approach has been selected. 901 902 Input Parameter: 903 + pc - the preconditioner contex 904 905 Application Interface Routine: PCPreSolve() 906 907 Notes: 908 The interface routine PCPreSolve() is not usually called directly by 909 the user, but instead is called by KSPSolve(). 910 */ 911 static PetscErrorCode PCPreSolve_BDDC(PC pc, KSP ksp, Vec rhs, Vec x) 912 { 913 PetscErrorCode ierr; 914 PC_BDDC *pcbddc = (PC_BDDC*)pc->data; 915 PC_IS *pcis = (PC_IS*)(pc->data); 916 Mat_IS *matis = (Mat_IS*)pc->pmat->data; 917 Mat temp_mat; 918 IS dirIS; 919 Vec used_vec; 920 PetscBool guess_nonzero; 921 922 PetscFunctionBegin; 923 /* if we are working with cg, one dirichlet solve can be avoided during Krylov iterations */ 924 if (ksp) { 925 PetscBool iscg; 926 ierr = PetscObjectTypeCompare((PetscObject)ksp,KSPCG,&iscg);CHKERRQ(ierr); 927 if (!iscg) { 928 ierr = PCBDDCSetUseExactDirichlet(pc,PETSC_FALSE);CHKERRQ(ierr); 929 } 930 } 931 /* Creates parallel work vectors used in presolve */ 932 if (!pcbddc->original_rhs) { 933 ierr = VecDuplicate(pcis->vec1_global,&pcbddc->original_rhs);CHKERRQ(ierr); 934 } 935 if (!pcbddc->temp_solution) { 936 ierr = VecDuplicate(pcis->vec1_global,&pcbddc->temp_solution);CHKERRQ(ierr); 937 } 938 if (x) { 939 ierr = PetscObjectReference((PetscObject)x);CHKERRQ(ierr); 940 used_vec = x; 941 } else { 942 ierr = PetscObjectReference((PetscObject)pcbddc->temp_solution);CHKERRQ(ierr); 943 used_vec = pcbddc->temp_solution; 944 ierr = VecSet(used_vec,0.0);CHKERRQ(ierr); 945 } 946 /* hack into ksp data structure PCPreSolve comes earlier in src/ksp/ksp/interface/itfunc.c */ 947 if (ksp) { 948 ierr = KSPGetInitialGuessNonzero(ksp,&guess_nonzero);CHKERRQ(ierr); 949 if (!guess_nonzero) { 950 ierr = VecSet(used_vec,0.0);CHKERRQ(ierr); 951 } 952 } 953 954 /* store the original rhs */ 955 ierr = VecCopy(rhs,pcbddc->original_rhs);CHKERRQ(ierr); 956 957 /* Take into account zeroed rows -> change rhs and store solution removed */ 958 /* note that Dirichlet boundaries in global ordering (if any) has already been translated into local ordering in PCBDDCAnalyzeInterface */ 959 ierr = PCBDDCGetDirichletBoundariesLocal(pc,&dirIS);CHKERRQ(ierr); 960 if (rhs && dirIS) { 961 PetscInt dirsize,i,*is_indices; 962 PetscScalar *array_x,*array_diagonal; 963 964 ierr = MatGetDiagonal(pc->pmat,pcis->vec1_global);CHKERRQ(ierr); 965 ierr = VecPointwiseDivide(pcis->vec1_global,rhs,pcis->vec1_global);CHKERRQ(ierr); 966 ierr = VecScatterBegin(matis->ctx,pcis->vec1_global,pcis->vec2_N,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 967 ierr = VecScatterEnd(matis->ctx,pcis->vec1_global,pcis->vec2_N,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 968 ierr = VecScatterBegin(matis->ctx,used_vec,pcis->vec1_N,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 969 ierr = VecScatterEnd(matis->ctx,used_vec,pcis->vec1_N,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 970 ierr = ISGetLocalSize(dirIS,&dirsize);CHKERRQ(ierr); 971 ierr = VecGetArray(pcis->vec1_N,&array_x);CHKERRQ(ierr); 972 ierr = VecGetArray(pcis->vec2_N,&array_diagonal);CHKERRQ(ierr); 973 ierr = ISGetIndices(dirIS,(const PetscInt**)&is_indices);CHKERRQ(ierr); 974 for (i=0; i<dirsize; i++) array_x[is_indices[i]] = array_diagonal[is_indices[i]]; 975 ierr = ISRestoreIndices(dirIS,(const PetscInt**)&is_indices);CHKERRQ(ierr); 976 ierr = VecRestoreArray(pcis->vec2_N,&array_diagonal);CHKERRQ(ierr); 977 ierr = VecRestoreArray(pcis->vec1_N,&array_x);CHKERRQ(ierr); 978 ierr = VecScatterBegin(matis->ctx,pcis->vec1_N,used_vec,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 979 ierr = VecScatterEnd(matis->ctx,pcis->vec1_N,used_vec,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 980 981 /* remove the computed solution from the rhs */ 982 ierr = VecScale(used_vec,-1.0);CHKERRQ(ierr); 983 ierr = MatMultAdd(pc->pmat,used_vec,rhs,rhs);CHKERRQ(ierr); 984 ierr = VecScale(used_vec,-1.0);CHKERRQ(ierr); 985 } 986 987 /* store partially computed solution and set initial guess */ 988 if (x) { 989 ierr = VecCopy(used_vec,pcbddc->temp_solution);CHKERRQ(ierr); 990 ierr = VecSet(used_vec,0.0);CHKERRQ(ierr); 991 if (pcbddc->use_exact_dirichlet_trick) { 992 ierr = VecScatterBegin(pcis->global_to_D,rhs,pcis->vec1_D,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 993 ierr = VecScatterEnd(pcis->global_to_D,rhs,pcis->vec1_D,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 994 ierr = KSPSolve(pcbddc->ksp_D,pcis->vec1_D,pcis->vec2_D);CHKERRQ(ierr); 995 ierr = VecScatterBegin(pcis->global_to_D,pcis->vec2_D,used_vec,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 996 ierr = VecScatterEnd(pcis->global_to_D,pcis->vec2_D,used_vec,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 997 if (ksp) { 998 ierr = KSPSetInitialGuessNonzero(ksp,PETSC_TRUE);CHKERRQ(ierr); 999 } 1000 } 1001 } 1002 1003 /* prepare MatMult and rhs for solver */ 1004 if (pcbddc->ChangeOfBasisMatrix) { 1005 /* swap pointers for local matrices */ 1006 temp_mat = matis->A; 1007 matis->A = pcbddc->local_mat; 1008 pcbddc->local_mat = temp_mat; 1009 if (rhs) { /* TODO the change of basis is fragile when the change involves dofs which are not on the same subset of the interface (i.e. Toselli's for curl-curl) */ 1010 /* Get local rhs and apply transformation of basis */ 1011 ierr = VecScatterBegin(pcis->global_to_B,rhs,pcis->vec1_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1012 ierr = VecScatterEnd(pcis->global_to_B,rhs,pcis->vec1_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1013 ierr = VecPointwiseMult(pcis->vec2_B,pcis->vec2_B,pcis->D);CHKERRQ(ierr); 1014 /* from original basis to modified basis */ 1015 ierr = MatMultTranspose(pcbddc->ChangeOfBasisMatrix,pcis->vec1_B,pcis->vec2_B);CHKERRQ(ierr); 1016 /* put back modified values into the global vec using ADD_VALUES copy mode */ 1017 ierr = VecScatterBegin(pcis->global_to_B,pcis->vec2_B,rhs,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1018 ierr = VecScatterEnd(pcis->global_to_B,pcis->vec2_B,rhs,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1019 } 1020 } 1021 1022 /* remove nullspace if present */ 1023 if (ksp && pcbddc->NullSpace) { 1024 ierr = MatNullSpaceRemove(pcbddc->NullSpace,used_vec);CHKERRQ(ierr); 1025 ierr = MatNullSpaceRemove(pcbddc->NullSpace,rhs);CHKERRQ(ierr); 1026 } 1027 ierr = VecDestroy(&used_vec);CHKERRQ(ierr); 1028 PetscFunctionReturn(0); 1029 } 1030 /* -------------------------------------------------------------------------- */ 1031 #undef __FUNCT__ 1032 #define __FUNCT__ "PCPostSolve_BDDC" 1033 /* -------------------------------------------------------------------------- */ 1034 /* 1035 PCPostSolve_BDDC - Changes the computed solution if a transformation of basis 1036 approach has been selected. Also, restores rhs to its original state. 1037 1038 Input Parameter: 1039 + pc - the preconditioner contex 1040 1041 Application Interface Routine: PCPostSolve() 1042 1043 Notes: 1044 The interface routine PCPostSolve() is not usually called directly by 1045 the user, but instead is called by KSPSolve(). 1046 */ 1047 static PetscErrorCode PCPostSolve_BDDC(PC pc, KSP ksp, Vec rhs, Vec x) 1048 { 1049 PetscErrorCode ierr; 1050 PC_BDDC *pcbddc = (PC_BDDC*)pc->data; 1051 PC_IS *pcis = (PC_IS*)(pc->data); 1052 Mat_IS *matis = (Mat_IS*)pc->pmat->data; 1053 Mat temp_mat; 1054 1055 PetscFunctionBegin; 1056 if (pcbddc->ChangeOfBasisMatrix) { 1057 /* swap pointers for local matrices */ 1058 temp_mat = matis->A; 1059 matis->A = pcbddc->local_mat; 1060 pcbddc->local_mat = temp_mat; 1061 } 1062 if (pcbddc->ChangeOfBasisMatrix && x) { 1063 /* Get Local boundary and apply transformation of basis to solution vector */ 1064 ierr = VecScatterBegin(pcis->global_to_B,x,pcis->vec1_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1065 ierr = VecScatterEnd(pcis->global_to_B,x,pcis->vec1_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1066 /* from modified basis to original basis */ 1067 ierr = MatMult(pcbddc->ChangeOfBasisMatrix,pcis->vec1_B,pcis->vec2_B);CHKERRQ(ierr); 1068 /* put back modified values into the global vec using INSERT_VALUES copy mode */ 1069 ierr = VecScatterBegin(pcis->global_to_B,pcis->vec2_B,x,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1070 ierr = VecScatterEnd(pcis->global_to_B,pcis->vec2_B,x,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1071 } 1072 /* add solution removed in presolve */ 1073 if (x) { 1074 ierr = VecAXPY(x,1.0,pcbddc->temp_solution);CHKERRQ(ierr); 1075 } 1076 /* restore rhs to its original state */ 1077 if (rhs) { 1078 ierr = VecCopy(pcbddc->original_rhs,rhs);CHKERRQ(ierr); 1079 } 1080 PetscFunctionReturn(0); 1081 } 1082 /* -------------------------------------------------------------------------- */ 1083 #undef __FUNCT__ 1084 #define __FUNCT__ "PCSetUp_BDDC" 1085 /* -------------------------------------------------------------------------- */ 1086 /* 1087 PCSetUp_BDDC - Prepares for the use of the BDDC preconditioner 1088 by setting data structures and options. 1089 1090 Input Parameter: 1091 + pc - the preconditioner context 1092 1093 Application Interface Routine: PCSetUp() 1094 1095 Notes: 1096 The interface routine PCSetUp() is not usually called directly by 1097 the user, but instead is called by PCApply() if necessary. 1098 */ 1099 PetscErrorCode PCSetUp_BDDC(PC pc) 1100 { 1101 PetscErrorCode ierr; 1102 PC_BDDC* pcbddc = (PC_BDDC*)pc->data; 1103 MatNullSpace nearnullspace; 1104 PetscBool computeis,computetopography,computesolvers; 1105 PetscBool new_nearnullspace_provided; 1106 1107 PetscFunctionBegin; 1108 /* the following lines of code should be replaced by a better logic between PCIS, PCNN, PCBDDC and other future nonoverlapping preconditioners */ 1109 /* For BDDC we need to define a local "Neumann" problem different to that defined in PCISSetup 1110 Also, BDDC directly build the Dirichlet problem */ 1111 /* split work */ 1112 if (pc->setupcalled) { 1113 computeis = PETSC_FALSE; 1114 if (pc->flag == SAME_NONZERO_PATTERN) { 1115 computetopography = PETSC_FALSE; 1116 computesolvers = PETSC_TRUE; 1117 } else { /* DIFFERENT_NONZERO_PATTERN */ 1118 computetopography = PETSC_TRUE; 1119 computesolvers = PETSC_TRUE; 1120 } 1121 } else { 1122 computeis = PETSC_TRUE; 1123 computetopography = PETSC_TRUE; 1124 computesolvers = PETSC_TRUE; 1125 } 1126 if (pcbddc->recompute_topography) { 1127 computetopography = PETSC_TRUE; 1128 } 1129 1130 /* Get stdout for dbg */ 1131 if (pcbddc->dbg_flag) { 1132 if (!pcbddc->dbg_viewer) { 1133 pcbddc->dbg_viewer = PETSC_VIEWER_STDOUT_(PetscObjectComm((PetscObject)pc)); 1134 ierr = PetscViewerASCIISynchronizedAllow(pcbddc->dbg_viewer,PETSC_TRUE);CHKERRQ(ierr); 1135 } 1136 ierr = PetscViewerASCIIAddTab(pcbddc->dbg_viewer,2*pcbddc->current_level);CHKERRQ(ierr); 1137 } 1138 1139 /* Set up all the "iterative substructuring" common block without computing solvers */ 1140 if (computeis) { 1141 /* HACK INTO PCIS */ 1142 PC_IS* pcis = (PC_IS*)pc->data; 1143 pcis->computesolvers = PETSC_FALSE; 1144 ierr = PCISSetUp(pc);CHKERRQ(ierr); 1145 ierr = ISLocalToGlobalMappingCreateIS(pcis->is_B_local,&pcbddc->BtoNmap);CHKERRQ(ierr); 1146 } 1147 1148 /* check user defined change of basis (if any) */ 1149 if (pcbddc->user_ChangeOfBasisMatrix) { 1150 PC_IS* pcis= (PC_IS*)pc->data; 1151 PetscInt n,m; 1152 ierr = MatGetSize(pcbddc->user_ChangeOfBasisMatrix,&n,&m);CHKERRQ(ierr); 1153 if (n != m) { 1154 SETERRQ2(PETSC_COMM_WORLD,PETSC_ERR_SUP,"Change of basis matrix should be square %d != %d\n",n,m); 1155 } else if (n != pcis->n_B && n != pcis->n) { 1156 SETERRQ3(PETSC_COMM_WORLD,PETSC_ERR_SUP,"Size of change of basis matrix %d differs from boundary size %d and local size %d\n",n,pcis->n_B,pcis->n); 1157 } 1158 /* use_change_of_basis flag is used to automatically compute a change of basis from constraints */ 1159 pcbddc->use_change_of_basis = PETSC_FALSE; 1160 } 1161 /* Analyze interface */ 1162 if (computetopography) { 1163 ierr = PCBDDCAnalyzeInterface(pc);CHKERRQ(ierr); 1164 /* Schurs on subsets should be reset */ 1165 if (pcbddc->deluxe_ctx) { 1166 ierr = PCBDDCSubSchursReset(pcbddc->deluxe_ctx->sub_schurs);CHKERRQ(ierr); 1167 } 1168 } 1169 1170 /* infer if NullSpace object attached to Mat via MatSetNearNullSpace has changed */ 1171 new_nearnullspace_provided = PETSC_FALSE; 1172 ierr = MatGetNearNullSpace(pc->pmat,&nearnullspace);CHKERRQ(ierr); 1173 if (pcbddc->onearnullspace) { /* already used nearnullspace */ 1174 if (!nearnullspace) { /* near null space attached to mat has been destroyed */ 1175 new_nearnullspace_provided = PETSC_TRUE; 1176 } else { 1177 /* determine if the two nullspaces are different (should be lightweight) */ 1178 if (nearnullspace != pcbddc->onearnullspace) { 1179 new_nearnullspace_provided = PETSC_TRUE; 1180 } else { /* maybe the user has changed the content of the nearnullspace so check vectors ObjectStateId */ 1181 PetscInt i; 1182 const Vec *nearnullvecs; 1183 PetscObjectState state; 1184 PetscInt nnsp_size; 1185 ierr = MatNullSpaceGetVecs(nearnullspace,NULL,&nnsp_size,&nearnullvecs);CHKERRQ(ierr); 1186 for (i=0;i<nnsp_size;i++) { 1187 ierr = PetscObjectStateGet((PetscObject)nearnullvecs[i],&state);CHKERRQ(ierr); 1188 if (pcbddc->onearnullvecs_state[i] != state) { 1189 new_nearnullspace_provided = PETSC_TRUE; 1190 break; 1191 } 1192 } 1193 } 1194 } 1195 } else { 1196 if (!nearnullspace) { /* both nearnullspaces are null */ 1197 new_nearnullspace_provided = PETSC_FALSE; 1198 } else { /* nearnullspace attached later */ 1199 new_nearnullspace_provided = PETSC_TRUE; 1200 } 1201 } 1202 1203 /* Setup constraints and related work vectors */ 1204 /* reset primal space flags */ 1205 pcbddc->new_primal_space = PETSC_FALSE; 1206 pcbddc->new_primal_space_local = PETSC_FALSE; 1207 if (computetopography || new_nearnullspace_provided) { 1208 /* It also sets the primal space flags */ 1209 ierr = PCBDDCConstraintsSetUp(pc);CHKERRQ(ierr); 1210 /* Allocate needed local vectors (which depends on quantities defined during ConstraintsSetUp) */ 1211 ierr = PCBDDCSetUpLocalWorkVectors(pc);CHKERRQ(ierr); 1212 /* Schurs on subsets should be reset */ 1213 if (pcbddc->deluxe_ctx) { 1214 ierr = PCBDDCSubSchursReset(pcbddc->deluxe_ctx->sub_schurs);CHKERRQ(ierr); 1215 } 1216 } 1217 1218 if (computesolvers || pcbddc->new_primal_space) { 1219 /* Create coarse and local stuffs */ 1220 ierr = PCBDDCSetUpSolvers(pc);CHKERRQ(ierr); 1221 /* Create scaling operators */ 1222 ierr = PCBDDCScalingSetUp(pc);CHKERRQ(ierr); 1223 } 1224 1225 if (pcbddc->dbg_flag) { 1226 ierr = PetscViewerASCIISubtractTab(pcbddc->dbg_viewer,2*pcbddc->current_level);CHKERRQ(ierr); 1227 } 1228 PetscFunctionReturn(0); 1229 } 1230 1231 /* -------------------------------------------------------------------------- */ 1232 /* 1233 PCApply_BDDC - Applies the BDDC operator to a vector. 1234 1235 Input Parameters: 1236 . pc - the preconditioner context 1237 . r - input vector (global) 1238 1239 Output Parameter: 1240 . z - output vector (global) 1241 1242 Application Interface Routine: PCApply() 1243 */ 1244 #undef __FUNCT__ 1245 #define __FUNCT__ "PCApply_BDDC" 1246 PetscErrorCode PCApply_BDDC(PC pc,Vec r,Vec z) 1247 { 1248 PC_IS *pcis = (PC_IS*)(pc->data); 1249 PC_BDDC *pcbddc = (PC_BDDC*)(pc->data); 1250 PetscErrorCode ierr; 1251 const PetscScalar one = 1.0; 1252 const PetscScalar m_one = -1.0; 1253 const PetscScalar zero = 0.0; 1254 1255 /* This code is similar to that provided in nn.c for PCNN 1256 NN interface preconditioner changed to BDDC 1257 Added support for M_3 preconditioner in the reference article (code is active if pcbddc->switch_static = PETSC_TRUE) */ 1258 1259 PetscFunctionBegin; 1260 if (!pcbddc->use_exact_dirichlet_trick) { 1261 /* First Dirichlet solve */ 1262 ierr = VecScatterBegin(pcis->global_to_D,r,pcis->vec1_D,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1263 ierr = VecScatterEnd(pcis->global_to_D,r,pcis->vec1_D,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1264 ierr = KSPSolve(pcbddc->ksp_D,pcis->vec1_D,pcis->vec2_D);CHKERRQ(ierr); 1265 /* 1266 Assembling right hand side for BDDC operator 1267 - pcis->vec1_D for the Dirichlet part (if needed, i.e. prec_flag=PETSC_TRUE) 1268 - pcis->vec1_B the interface part of the global vector z 1269 */ 1270 ierr = VecScale(pcis->vec2_D,m_one);CHKERRQ(ierr); 1271 ierr = MatMult(pcis->A_BI,pcis->vec2_D,pcis->vec1_B);CHKERRQ(ierr); 1272 if (pcbddc->switch_static) { ierr = MatMultAdd(pcis->A_II,pcis->vec2_D,pcis->vec1_D,pcis->vec1_D);CHKERRQ(ierr); } 1273 ierr = VecScale(pcis->vec2_D,m_one);CHKERRQ(ierr); 1274 ierr = VecCopy(r,z);CHKERRQ(ierr); 1275 ierr = VecScatterBegin(pcis->global_to_B,pcis->vec1_B,z,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1276 ierr = VecScatterEnd(pcis->global_to_B,pcis->vec1_B,z,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1277 ierr = PCBDDCScalingRestriction(pc,z,pcis->vec1_B);CHKERRQ(ierr); 1278 } else { 1279 ierr = VecSet(pcis->vec1_D,zero);CHKERRQ(ierr); 1280 ierr = VecSet(pcis->vec2_D,zero);CHKERRQ(ierr); 1281 ierr = PCBDDCScalingRestriction(pc,r,pcis->vec1_B);CHKERRQ(ierr); 1282 } 1283 1284 /* Apply interface preconditioner 1285 input/output vecs: pcis->vec1_B and pcis->vec1_D */ 1286 ierr = PCBDDCApplyInterfacePreconditioner(pc,PETSC_FALSE);CHKERRQ(ierr); 1287 1288 /* Apply transpose of partition of unity operator */ 1289 ierr = PCBDDCScalingExtension(pc,pcis->vec1_B,z);CHKERRQ(ierr); 1290 1291 /* Second Dirichlet solve and assembling of output */ 1292 ierr = VecScatterBegin(pcis->global_to_B,z,pcis->vec1_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1293 ierr = VecScatterEnd(pcis->global_to_B,z,pcis->vec1_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1294 ierr = MatMult(pcis->A_IB,pcis->vec1_B,pcis->vec3_D);CHKERRQ(ierr); 1295 if (pcbddc->switch_static) { ierr = MatMultAdd(pcis->A_II,pcis->vec1_D,pcis->vec3_D,pcis->vec3_D);CHKERRQ(ierr); } 1296 ierr = KSPSolve(pcbddc->ksp_D,pcis->vec3_D,pcis->vec4_D);CHKERRQ(ierr); 1297 ierr = VecScale(pcis->vec4_D,m_one);CHKERRQ(ierr); 1298 if (pcbddc->switch_static) { ierr = VecAXPY(pcis->vec4_D,one,pcis->vec1_D);CHKERRQ(ierr); } 1299 ierr = VecAXPY(pcis->vec2_D,one,pcis->vec4_D);CHKERRQ(ierr); 1300 ierr = VecScatterBegin(pcis->global_to_D,pcis->vec2_D,z,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1301 ierr = VecScatterEnd(pcis->global_to_D,pcis->vec2_D,z,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1302 PetscFunctionReturn(0); 1303 } 1304 1305 /* -------------------------------------------------------------------------- */ 1306 /* 1307 PCApplyTranspose_BDDC - Applies the transpose of the BDDC operator to a vector. 1308 1309 Input Parameters: 1310 . pc - the preconditioner context 1311 . r - input vector (global) 1312 1313 Output Parameter: 1314 . z - output vector (global) 1315 1316 Application Interface Routine: PCApplyTranspose() 1317 */ 1318 #undef __FUNCT__ 1319 #define __FUNCT__ "PCApplyTranspose_BDDC" 1320 PetscErrorCode PCApplyTranspose_BDDC(PC pc,Vec r,Vec z) 1321 { 1322 PC_IS *pcis = (PC_IS*)(pc->data); 1323 PC_BDDC *pcbddc = (PC_BDDC*)(pc->data); 1324 PetscErrorCode ierr; 1325 const PetscScalar one = 1.0; 1326 const PetscScalar m_one = -1.0; 1327 const PetscScalar zero = 0.0; 1328 1329 PetscFunctionBegin; 1330 if (!pcbddc->use_exact_dirichlet_trick) { 1331 /* First Dirichlet solve */ 1332 ierr = VecScatterBegin(pcis->global_to_D,r,pcis->vec1_D,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1333 ierr = VecScatterEnd(pcis->global_to_D,r,pcis->vec1_D,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1334 ierr = KSPSolveTranspose(pcbddc->ksp_D,pcis->vec1_D,pcis->vec2_D);CHKERRQ(ierr); 1335 /* 1336 Assembling right hand side for BDDC operator 1337 - pcis->vec1_D for the Dirichlet part (if needed, i.e. prec_flag=PETSC_TRUE) 1338 - pcis->vec1_B the interface part of the global vector z 1339 */ 1340 ierr = VecScale(pcis->vec2_D,m_one);CHKERRQ(ierr); 1341 ierr = MatMultTranspose(pcis->A_IB,pcis->vec2_D,pcis->vec1_B);CHKERRQ(ierr); 1342 if (pcbddc->switch_static) { ierr = MatMultTransposeAdd(pcis->A_II,pcis->vec2_D,pcis->vec1_D,pcis->vec1_D);CHKERRQ(ierr); } 1343 ierr = VecScale(pcis->vec2_D,m_one);CHKERRQ(ierr); 1344 ierr = VecCopy(r,z);CHKERRQ(ierr); 1345 ierr = VecScatterBegin(pcis->global_to_B,pcis->vec1_B,z,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1346 ierr = VecScatterEnd(pcis->global_to_B,pcis->vec1_B,z,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1347 ierr = PCBDDCScalingRestriction(pc,z,pcis->vec1_B);CHKERRQ(ierr); 1348 } else { 1349 ierr = VecSet(pcis->vec1_D,zero);CHKERRQ(ierr); 1350 ierr = VecSet(pcis->vec2_D,zero);CHKERRQ(ierr); 1351 ierr = PCBDDCScalingRestriction(pc,r,pcis->vec1_B);CHKERRQ(ierr); 1352 } 1353 1354 /* Apply interface preconditioner 1355 input/output vecs: pcis->vec1_B and pcis->vec1_D */ 1356 ierr = PCBDDCApplyInterfacePreconditioner(pc,PETSC_TRUE);CHKERRQ(ierr); 1357 1358 /* Apply transpose of partition of unity operator */ 1359 ierr = PCBDDCScalingExtension(pc,pcis->vec1_B,z);CHKERRQ(ierr); 1360 1361 /* Second Dirichlet solve and assembling of output */ 1362 ierr = VecScatterBegin(pcis->global_to_B,z,pcis->vec1_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1363 ierr = VecScatterEnd(pcis->global_to_B,z,pcis->vec1_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1364 ierr = MatMultTranspose(pcis->A_BI,pcis->vec1_B,pcis->vec3_D);CHKERRQ(ierr); 1365 if (pcbddc->switch_static) { ierr = MatMultTransposeAdd(pcis->A_II,pcis->vec1_D,pcis->vec3_D,pcis->vec3_D);CHKERRQ(ierr); } 1366 ierr = KSPSolveTranspose(pcbddc->ksp_D,pcis->vec3_D,pcis->vec4_D);CHKERRQ(ierr); 1367 ierr = VecScale(pcis->vec4_D,m_one);CHKERRQ(ierr); 1368 if (pcbddc->switch_static) { ierr = VecAXPY(pcis->vec4_D,one,pcis->vec1_D);CHKERRQ(ierr); } 1369 ierr = VecAXPY(pcis->vec2_D,one,pcis->vec4_D);CHKERRQ(ierr); 1370 ierr = VecScatterBegin(pcis->global_to_D,pcis->vec2_D,z,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1371 ierr = VecScatterEnd(pcis->global_to_D,pcis->vec2_D,z,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1372 PetscFunctionReturn(0); 1373 } 1374 /* -------------------------------------------------------------------------- */ 1375 1376 #undef __FUNCT__ 1377 #define __FUNCT__ "PCDestroy_BDDC" 1378 PetscErrorCode PCDestroy_BDDC(PC pc) 1379 { 1380 PC_BDDC *pcbddc = (PC_BDDC*)pc->data; 1381 PetscErrorCode ierr; 1382 1383 PetscFunctionBegin; 1384 /* free data created by PCIS */ 1385 ierr = PCISDestroy(pc);CHKERRQ(ierr); 1386 /* free BDDC custom data */ 1387 ierr = PCBDDCResetCustomization(pc);CHKERRQ(ierr); 1388 /* destroy objects related to topography */ 1389 ierr = PCBDDCResetTopography(pc);CHKERRQ(ierr); 1390 /* free allocated graph structure */ 1391 ierr = PetscFree(pcbddc->mat_graph);CHKERRQ(ierr); 1392 /* destroy objects for scaling operator */ 1393 ierr = PCBDDCScalingDestroy(pc);CHKERRQ(ierr); 1394 ierr = PetscFree(pcbddc->deluxe_ctx);CHKERRQ(ierr); 1395 /* free solvers stuff */ 1396 ierr = PCBDDCResetSolvers(pc);CHKERRQ(ierr); 1397 /* free global vectors needed in presolve */ 1398 ierr = VecDestroy(&pcbddc->temp_solution);CHKERRQ(ierr); 1399 ierr = VecDestroy(&pcbddc->original_rhs);CHKERRQ(ierr); 1400 ierr = ISLocalToGlobalMappingDestroy(&pcbddc->BtoNmap);CHKERRQ(ierr); 1401 /* remove functions */ 1402 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetChangeOfBasisLocalMat_C",NULL);CHKERRQ(ierr); 1403 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetPrimalVerticesLocalIS_C",NULL);CHKERRQ(ierr); 1404 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetCoarseningRatio_C",NULL);CHKERRQ(ierr); 1405 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetLevel_C",NULL);CHKERRQ(ierr); 1406 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetUseExactDirichlet_C",NULL);CHKERRQ(ierr); 1407 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetLevels_C",NULL);CHKERRQ(ierr); 1408 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetNullSpace_C",NULL);CHKERRQ(ierr); 1409 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetDirichletBoundaries_C",NULL);CHKERRQ(ierr); 1410 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetDirichletBoundariesLocal_C",NULL);CHKERRQ(ierr); 1411 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetNeumannBoundaries_C",NULL);CHKERRQ(ierr); 1412 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetNeumannBoundariesLocal_C",NULL);CHKERRQ(ierr); 1413 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCGetDirichletBoundaries_C",NULL);CHKERRQ(ierr); 1414 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCGetDirichletBoundariesLocal_C",NULL);CHKERRQ(ierr); 1415 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCGetNeumannBoundaries_C",NULL);CHKERRQ(ierr); 1416 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCGetNeumannBoundariesLocal_C",NULL);CHKERRQ(ierr); 1417 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetDofsSplitting_C",NULL);CHKERRQ(ierr); 1418 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetDofsSplittingLocal_C",NULL);CHKERRQ(ierr); 1419 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetLocalAdjacencyGraph_C",NULL);CHKERRQ(ierr); 1420 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCCreateFETIDPOperators_C",NULL);CHKERRQ(ierr); 1421 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCMatFETIDPGetRHS_C",NULL);CHKERRQ(ierr); 1422 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCMatFETIDPGetSolution_C",NULL);CHKERRQ(ierr); 1423 /* Free the private data structure */ 1424 ierr = PetscFree(pc->data);CHKERRQ(ierr); 1425 PetscFunctionReturn(0); 1426 } 1427 /* -------------------------------------------------------------------------- */ 1428 1429 #undef __FUNCT__ 1430 #define __FUNCT__ "PCBDDCMatFETIDPGetRHS_BDDC" 1431 static PetscErrorCode PCBDDCMatFETIDPGetRHS_BDDC(Mat fetidp_mat, Vec standard_rhs, Vec fetidp_flux_rhs) 1432 { 1433 FETIDPMat_ctx mat_ctx; 1434 PC_IS* pcis; 1435 PC_BDDC* pcbddc; 1436 PetscErrorCode ierr; 1437 1438 PetscFunctionBegin; 1439 ierr = MatShellGetContext(fetidp_mat,&mat_ctx);CHKERRQ(ierr); 1440 pcis = (PC_IS*)mat_ctx->pc->data; 1441 pcbddc = (PC_BDDC*)mat_ctx->pc->data; 1442 1443 /* change of basis for physical rhs if needed 1444 It also changes the rhs in case of dirichlet boundaries */ 1445 ierr = PCPreSolve_BDDC(mat_ctx->pc,NULL,standard_rhs,NULL);CHKERRQ(ierr); 1446 /* store vectors for computation of fetidp final solution */ 1447 ierr = VecScatterBegin(pcis->global_to_D,standard_rhs,mat_ctx->temp_solution_D,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1448 ierr = VecScatterEnd(pcis->global_to_D,standard_rhs,mat_ctx->temp_solution_D,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1449 /* scale rhs since it should be unassembled */ 1450 /* TODO use counter scaling? (also below) */ 1451 ierr = VecScatterBegin(pcis->global_to_B,standard_rhs,mat_ctx->temp_solution_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1452 ierr = VecScatterEnd(pcis->global_to_B,standard_rhs,mat_ctx->temp_solution_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1453 /* Apply partition of unity */ 1454 ierr = VecPointwiseMult(mat_ctx->temp_solution_B,pcis->D,mat_ctx->temp_solution_B);CHKERRQ(ierr); 1455 /* ierr = PCBDDCScalingRestriction(mat_ctx->pc,standard_rhs,mat_ctx->temp_solution_B);CHKERRQ(ierr); */ 1456 if (!pcbddc->switch_static) { 1457 /* compute partially subassembled Schur complement right-hand side */ 1458 ierr = KSPSolve(pcbddc->ksp_D,mat_ctx->temp_solution_D,pcis->vec1_D);CHKERRQ(ierr); 1459 ierr = MatMult(pcis->A_BI,pcis->vec1_D,pcis->vec1_B);CHKERRQ(ierr); 1460 ierr = VecAXPY(mat_ctx->temp_solution_B,-1.0,pcis->vec1_B);CHKERRQ(ierr); 1461 ierr = VecSet(standard_rhs,0.0);CHKERRQ(ierr); 1462 ierr = VecScatterBegin(pcis->global_to_B,mat_ctx->temp_solution_B,standard_rhs,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1463 ierr = VecScatterEnd(pcis->global_to_B,mat_ctx->temp_solution_B,standard_rhs,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1464 /* ierr = PCBDDCScalingRestriction(mat_ctx->pc,standard_rhs,mat_ctx->temp_solution_B);CHKERRQ(ierr); */ 1465 ierr = VecScatterBegin(pcis->global_to_B,standard_rhs,mat_ctx->temp_solution_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1466 ierr = VecScatterEnd(pcis->global_to_B,standard_rhs,mat_ctx->temp_solution_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1467 ierr = VecPointwiseMult(mat_ctx->temp_solution_B,pcis->D,mat_ctx->temp_solution_B);CHKERRQ(ierr); 1468 } 1469 /* BDDC rhs */ 1470 ierr = VecCopy(mat_ctx->temp_solution_B,pcis->vec1_B);CHKERRQ(ierr); 1471 if (pcbddc->switch_static) { 1472 ierr = VecCopy(mat_ctx->temp_solution_D,pcis->vec1_D);CHKERRQ(ierr); 1473 } 1474 /* apply BDDC */ 1475 ierr = PCBDDCApplyInterfacePreconditioner(mat_ctx->pc,PETSC_FALSE);CHKERRQ(ierr); 1476 /* Application of B_delta and assembling of rhs for fetidp fluxes */ 1477 ierr = VecSet(fetidp_flux_rhs,0.0);CHKERRQ(ierr); 1478 ierr = MatMult(mat_ctx->B_delta,pcis->vec1_B,mat_ctx->lambda_local);CHKERRQ(ierr); 1479 ierr = VecScatterBegin(mat_ctx->l2g_lambda,mat_ctx->lambda_local,fetidp_flux_rhs,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1480 ierr = VecScatterEnd (mat_ctx->l2g_lambda,mat_ctx->lambda_local,fetidp_flux_rhs,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1481 /* restore original rhs */ 1482 ierr = VecCopy(pcbddc->original_rhs,standard_rhs);CHKERRQ(ierr); 1483 PetscFunctionReturn(0); 1484 } 1485 1486 #undef __FUNCT__ 1487 #define __FUNCT__ "PCBDDCMatFETIDPGetRHS" 1488 /*@ 1489 PCBDDCMatFETIDPGetRHS - Compute the right-hand side for FETIDP linear system 1490 1491 Collective 1492 1493 Input Parameters: 1494 + fetidp_mat - the FETIDP matrix object obtained by calling PCBDDCCreateFETIDPOperators 1495 . standard_rhs - the right-hand side for your linear system 1496 1497 Output Parameters: 1498 - fetidp_flux_rhs - the right-hand side for the FETIDP linear system 1499 1500 Level: developer 1501 1502 Notes: 1503 1504 .seealso: PCBDDC,PCBDDCCreateFETIDPOperators 1505 @*/ 1506 PetscErrorCode PCBDDCMatFETIDPGetRHS(Mat fetidp_mat, Vec standard_rhs, Vec fetidp_flux_rhs) 1507 { 1508 FETIDPMat_ctx mat_ctx; 1509 PetscErrorCode ierr; 1510 1511 PetscFunctionBegin; 1512 ierr = MatShellGetContext(fetidp_mat,&mat_ctx);CHKERRQ(ierr); 1513 ierr = PetscTryMethod(mat_ctx->pc,"PCBDDCMatFETIDPGetRHS_C",(Mat,Vec,Vec),(fetidp_mat,standard_rhs,fetidp_flux_rhs));CHKERRQ(ierr); 1514 PetscFunctionReturn(0); 1515 } 1516 /* -------------------------------------------------------------------------- */ 1517 1518 #undef __FUNCT__ 1519 #define __FUNCT__ "PCBDDCMatFETIDPGetSolution_BDDC" 1520 static PetscErrorCode PCBDDCMatFETIDPGetSolution_BDDC(Mat fetidp_mat, Vec fetidp_flux_sol, Vec standard_sol) 1521 { 1522 FETIDPMat_ctx mat_ctx; 1523 PC_IS* pcis; 1524 PC_BDDC* pcbddc; 1525 PetscErrorCode ierr; 1526 1527 PetscFunctionBegin; 1528 ierr = MatShellGetContext(fetidp_mat,&mat_ctx);CHKERRQ(ierr); 1529 pcis = (PC_IS*)mat_ctx->pc->data; 1530 pcbddc = (PC_BDDC*)mat_ctx->pc->data; 1531 1532 /* apply B_delta^T */ 1533 ierr = VecScatterBegin(mat_ctx->l2g_lambda,fetidp_flux_sol,mat_ctx->lambda_local,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1534 ierr = VecScatterEnd (mat_ctx->l2g_lambda,fetidp_flux_sol,mat_ctx->lambda_local,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1535 ierr = MatMultTranspose(mat_ctx->B_delta,mat_ctx->lambda_local,pcis->vec1_B);CHKERRQ(ierr); 1536 /* compute rhs for BDDC application */ 1537 ierr = VecAYPX(pcis->vec1_B,-1.0,mat_ctx->temp_solution_B);CHKERRQ(ierr); 1538 if (pcbddc->switch_static) { 1539 ierr = VecCopy(mat_ctx->temp_solution_D,pcis->vec1_D);CHKERRQ(ierr); 1540 } 1541 /* apply BDDC */ 1542 ierr = PCBDDCApplyInterfacePreconditioner(mat_ctx->pc,PETSC_FALSE);CHKERRQ(ierr); 1543 /* put values into standard global vector */ 1544 ierr = VecScatterBegin(pcis->global_to_B,pcis->vec1_B,standard_sol,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1545 ierr = VecScatterEnd (pcis->global_to_B,pcis->vec1_B,standard_sol,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1546 if (!pcbddc->switch_static) { 1547 /* compute values into the interior if solved for the partially subassembled Schur complement */ 1548 ierr = MatMult(pcis->A_IB,pcis->vec1_B,pcis->vec1_D);CHKERRQ(ierr); 1549 ierr = VecAXPY(mat_ctx->temp_solution_D,-1.0,pcis->vec1_D);CHKERRQ(ierr); 1550 ierr = KSPSolve(pcbddc->ksp_D,mat_ctx->temp_solution_D,pcis->vec1_D);CHKERRQ(ierr); 1551 } 1552 ierr = VecScatterBegin(pcis->global_to_D,pcis->vec1_D,standard_sol,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1553 ierr = VecScatterEnd (pcis->global_to_D,pcis->vec1_D,standard_sol,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1554 /* final change of basis if needed 1555 Is also sums the dirichlet part removed during RHS assembling */ 1556 ierr = PCPostSolve_BDDC(mat_ctx->pc,NULL,NULL,standard_sol);CHKERRQ(ierr); 1557 PetscFunctionReturn(0); 1558 } 1559 1560 #undef __FUNCT__ 1561 #define __FUNCT__ "PCBDDCMatFETIDPGetSolution" 1562 /*@ 1563 PCBDDCMatFETIDPGetSolution - Compute the physical solution from the solution of the FETIDP linear system 1564 1565 Collective 1566 1567 Input Parameters: 1568 + fetidp_mat - the FETIDP matrix obtained by calling PCBDDCCreateFETIDPOperators 1569 . fetidp_flux_sol - the solution of the FETIDP linear system 1570 1571 Output Parameters: 1572 - standard_sol - the solution defined on the physical domain 1573 1574 Level: developer 1575 1576 Notes: 1577 1578 .seealso: PCBDDC,PCBDDCCreateFETIDPOperators 1579 @*/ 1580 PetscErrorCode PCBDDCMatFETIDPGetSolution(Mat fetidp_mat, Vec fetidp_flux_sol, Vec standard_sol) 1581 { 1582 FETIDPMat_ctx mat_ctx; 1583 PetscErrorCode ierr; 1584 1585 PetscFunctionBegin; 1586 ierr = MatShellGetContext(fetidp_mat,&mat_ctx);CHKERRQ(ierr); 1587 ierr = PetscTryMethod(mat_ctx->pc,"PCBDDCMatFETIDPGetSolution_C",(Mat,Vec,Vec),(fetidp_mat,fetidp_flux_sol,standard_sol));CHKERRQ(ierr); 1588 PetscFunctionReturn(0); 1589 } 1590 /* -------------------------------------------------------------------------- */ 1591 1592 extern PetscErrorCode FETIDPMatMult(Mat,Vec,Vec); 1593 extern PetscErrorCode FETIDPMatMultTranspose(Mat,Vec,Vec); 1594 extern PetscErrorCode PCBDDCDestroyFETIDPMat(Mat); 1595 extern PetscErrorCode FETIDPPCApply(PC,Vec,Vec); 1596 extern PetscErrorCode FETIDPPCApplyTranspose(PC,Vec,Vec); 1597 extern PetscErrorCode PCBDDCDestroyFETIDPPC(PC); 1598 1599 #undef __FUNCT__ 1600 #define __FUNCT__ "PCBDDCCreateFETIDPOperators_BDDC" 1601 static PetscErrorCode PCBDDCCreateFETIDPOperators_BDDC(PC pc, Mat *fetidp_mat, PC *fetidp_pc) 1602 { 1603 1604 FETIDPMat_ctx fetidpmat_ctx; 1605 Mat newmat; 1606 FETIDPPC_ctx fetidppc_ctx; 1607 PC newpc; 1608 MPI_Comm comm; 1609 PetscErrorCode ierr; 1610 1611 PetscFunctionBegin; 1612 ierr = PetscObjectGetComm((PetscObject)pc,&comm);CHKERRQ(ierr); 1613 /* FETIDP linear matrix */ 1614 ierr = PCBDDCCreateFETIDPMatContext(pc,&fetidpmat_ctx);CHKERRQ(ierr); 1615 ierr = PCBDDCSetupFETIDPMatContext(fetidpmat_ctx);CHKERRQ(ierr); 1616 ierr = MatCreateShell(comm,PETSC_DECIDE,PETSC_DECIDE,fetidpmat_ctx->n_lambda,fetidpmat_ctx->n_lambda,fetidpmat_ctx,&newmat);CHKERRQ(ierr); 1617 ierr = MatShellSetOperation(newmat,MATOP_MULT,(void (*)(void))FETIDPMatMult);CHKERRQ(ierr); 1618 ierr = MatShellSetOperation(newmat,MATOP_MULT_TRANSPOSE,(void (*)(void))FETIDPMatMultTranspose);CHKERRQ(ierr); 1619 ierr = MatShellSetOperation(newmat,MATOP_DESTROY,(void (*)(void))PCBDDCDestroyFETIDPMat);CHKERRQ(ierr); 1620 ierr = MatSetUp(newmat);CHKERRQ(ierr); 1621 /* FETIDP preconditioner */ 1622 ierr = PCBDDCCreateFETIDPPCContext(pc,&fetidppc_ctx);CHKERRQ(ierr); 1623 ierr = PCBDDCSetupFETIDPPCContext(newmat,fetidppc_ctx);CHKERRQ(ierr); 1624 ierr = PCCreate(comm,&newpc);CHKERRQ(ierr); 1625 ierr = PCSetType(newpc,PCSHELL);CHKERRQ(ierr); 1626 ierr = PCShellSetContext(newpc,fetidppc_ctx);CHKERRQ(ierr); 1627 ierr = PCShellSetApply(newpc,FETIDPPCApply);CHKERRQ(ierr); 1628 ierr = PCShellSetApplyTranspose(newpc,FETIDPPCApplyTranspose);CHKERRQ(ierr); 1629 ierr = PCShellSetDestroy(newpc,PCBDDCDestroyFETIDPPC);CHKERRQ(ierr); 1630 ierr = PCSetOperators(newpc,newmat,newmat);CHKERRQ(ierr); 1631 ierr = PCSetUp(newpc);CHKERRQ(ierr); 1632 /* return pointers for objects created */ 1633 *fetidp_mat=newmat; 1634 *fetidp_pc=newpc; 1635 PetscFunctionReturn(0); 1636 } 1637 1638 #undef __FUNCT__ 1639 #define __FUNCT__ "PCBDDCCreateFETIDPOperators" 1640 /*@ 1641 PCBDDCCreateFETIDPOperators - Create operators for FETIDP 1642 1643 Collective 1644 1645 Input Parameters: 1646 + pc - the BDDC preconditioning context already setup 1647 1648 Output Parameters: 1649 . fetidp_mat - shell FETIDP matrix object 1650 . fetidp_pc - shell Dirichlet preconditioner for FETIDP matrix 1651 1652 Options Database Keys: 1653 - -fetidp_fullyredundant: use or not a fully redundant set of Lagrange multipliers 1654 1655 Level: developer 1656 1657 Notes: 1658 Currently the only operation provided for FETIDP matrix is MatMult 1659 1660 .seealso: PCBDDC 1661 @*/ 1662 PetscErrorCode PCBDDCCreateFETIDPOperators(PC pc, Mat *fetidp_mat, PC *fetidp_pc) 1663 { 1664 PetscErrorCode ierr; 1665 1666 PetscFunctionBegin; 1667 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1668 if (pc->setupcalled) { 1669 ierr = PetscUseMethod(pc,"PCBDDCCreateFETIDPOperators_C",(PC,Mat*,PC*),(pc,fetidp_mat,fetidp_pc));CHKERRQ(ierr); 1670 } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"You must call PCSetup_BDDC() first \n"); 1671 PetscFunctionReturn(0); 1672 } 1673 /* -------------------------------------------------------------------------- */ 1674 /*MC 1675 PCBDDC - Balancing Domain Decomposition by Constraints. 1676 1677 An implementation of the BDDC preconditioner based on 1678 1679 .vb 1680 [1] C. R. Dohrmann. "An approximate BDDC preconditioner", Numerical Linear Algebra with Applications Volume 14, Issue 2, pages 149-168, March 2007 1681 [2] A. Klawonn and O. B. Widlund. "Dual-Primal FETI Methods for Linear Elasticity", http://cs.nyu.edu/csweb/Research/TechReports/TR2004-855/TR2004-855.pdf 1682 [3] J. Mandel, B. Sousedik, C. R. Dohrmann. "Multispace and Multilevel BDDC", http://arxiv.org/abs/0712.3977 1683 .ve 1684 1685 The matrix to be preconditioned (Pmat) must be of type MATIS. 1686 1687 Currently works with MATIS matrices with local Neumann matrices of type MATSEQAIJ, MATSEQBAIJ or MATSEQSBAIJ, either with real or complex numbers. 1688 1689 It also works with unsymmetric and indefinite problems. 1690 1691 Unlike 'conventional' interface preconditioners, PCBDDC iterates over all degrees of freedom, not just those on the interface. This allows the use of approximate solvers on the subdomains. 1692 1693 Approximate local solvers are automatically adapted for singular linear problems (see [1]) if the user has provided the nullspace using PCBDDCSetNullSpace 1694 1695 Boundary nodes are split in vertices, edges and faces using information from the local to global mapping of dofs and the local connectivity graph of nodes. The latter can be customized by using PCBDDCSetLocalAdjacencyGraph() 1696 1697 Constraints can be customized by attaching a MatNullSpace object to the MATIS matrix via MatSetNearNullSpace(). 1698 1699 Change of basis is performed similarly to [2] when requested. When more the one constraint is present on a single connected component (i.e. an edge or a face), a robust method based on local QR factorizations is used. 1700 1701 The PETSc implementation also supports multilevel BDDC [3]. Coarse grids are partitioned using MatPartitioning object. 1702 1703 Options Database Keys: 1704 1705 . -pc_bddc_use_vertices <1> - use or not vertices in primal space 1706 . -pc_bddc_use_edges <1> - use or not edges in primal space 1707 . -pc_bddc_use_faces <0> - use or not faces in primal space 1708 . -pc_bddc_use_change_of_basis <0> - use change of basis approach (on edges only) 1709 . -pc_bddc_use_change_on_faces <0> - use change of basis approach on faces if change of basis has been requested 1710 . -pc_bddc_switch_static <0> - switches from M_2 to M_3 operator (see reference article [1]) 1711 . -pc_bddc_levels <0> - maximum number of levels for multilevel 1712 . -pc_bddc_coarsening_ratio - H/h ratio at the coarser level 1713 - -pc_bddc_check_level <0> - set verbosity level of debugging output 1714 1715 Options for Dirichlet, Neumann or coarse solver can be set with 1716 .vb 1717 -pc_bddc_dirichlet_ 1718 -pc_bddc_neumann_ 1719 -pc_bddc_coarse_ 1720 .ve 1721 e.g -pc_bddc_dirichlet_ksp_type richardson -pc_bddc_dirichlet_pc_type gamg 1722 1723 When using a multilevel approach, solvers' options at the N-th level can be specified as 1724 .vb 1725 -pc_bddc_dirichlet_lN_ 1726 -pc_bddc_neumann_lN_ 1727 -pc_bddc_coarse_lN_ 1728 .ve 1729 Note that level number ranges from the finest 0 to the coarsest N. 1730 1731 Level: intermediate 1732 1733 Developer notes: 1734 1735 New deluxe scaling operator will be available soon. 1736 1737 Contributed by Stefano Zampini 1738 1739 .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PC, MATIS 1740 M*/ 1741 1742 #undef __FUNCT__ 1743 #define __FUNCT__ "PCCreate_BDDC" 1744 PETSC_EXTERN PetscErrorCode PCCreate_BDDC(PC pc) 1745 { 1746 PetscErrorCode ierr; 1747 PC_BDDC *pcbddc; 1748 1749 PetscFunctionBegin; 1750 /* Creates the private data structure for this preconditioner and attach it to the PC object. */ 1751 ierr = PetscNewLog(pc,&pcbddc);CHKERRQ(ierr); 1752 pc->data = (void*)pcbddc; 1753 1754 /* create PCIS data structure */ 1755 ierr = PCISCreate(pc);CHKERRQ(ierr); 1756 1757 /* BDDC customization */ 1758 pcbddc->use_local_adj = PETSC_TRUE; 1759 pcbddc->use_vertices = PETSC_TRUE; 1760 pcbddc->use_edges = PETSC_TRUE; 1761 pcbddc->use_faces = PETSC_FALSE; 1762 pcbddc->use_change_of_basis = PETSC_FALSE; 1763 pcbddc->use_change_on_faces = PETSC_FALSE; 1764 pcbddc->switch_static = PETSC_FALSE; 1765 pcbddc->use_nnsp_true = PETSC_FALSE; /* not yet exposed */ 1766 pcbddc->dbg_flag = 0; 1767 /* private */ 1768 pcbddc->issym = PETSC_FALSE; 1769 pcbddc->BtoNmap = 0; 1770 pcbddc->local_primal_size = 0; 1771 pcbddc->n_vertices = 0; 1772 pcbddc->n_actual_vertices = 0; 1773 pcbddc->n_constraints = 0; 1774 pcbddc->primal_indices_local_idxs = 0; 1775 pcbddc->recompute_topography = PETSC_FALSE; 1776 pcbddc->coarse_size = -1; 1777 pcbddc->new_primal_space = PETSC_FALSE; 1778 pcbddc->new_primal_space_local = PETSC_FALSE; 1779 pcbddc->global_primal_indices = 0; 1780 pcbddc->onearnullspace = 0; 1781 pcbddc->onearnullvecs_state = 0; 1782 pcbddc->user_primal_vertices = 0; 1783 pcbddc->NullSpace = 0; 1784 pcbddc->temp_solution = 0; 1785 pcbddc->original_rhs = 0; 1786 pcbddc->local_mat = 0; 1787 pcbddc->ChangeOfBasisMatrix = 0; 1788 pcbddc->user_ChangeOfBasisMatrix = 0; 1789 pcbddc->coarse_vec = 0; 1790 pcbddc->coarse_rhs = 0; 1791 pcbddc->coarse_ksp = 0; 1792 pcbddc->coarse_phi_B = 0; 1793 pcbddc->coarse_phi_D = 0; 1794 pcbddc->coarse_psi_B = 0; 1795 pcbddc->coarse_psi_D = 0; 1796 pcbddc->vec1_P = 0; 1797 pcbddc->vec1_R = 0; 1798 pcbddc->vec2_R = 0; 1799 pcbddc->local_auxmat1 = 0; 1800 pcbddc->local_auxmat2 = 0; 1801 pcbddc->R_to_B = 0; 1802 pcbddc->R_to_D = 0; 1803 pcbddc->ksp_D = 0; 1804 pcbddc->ksp_R = 0; 1805 pcbddc->NeumannBoundaries = 0; 1806 pcbddc->NeumannBoundariesLocal = 0; 1807 pcbddc->DirichletBoundaries = 0; 1808 pcbddc->DirichletBoundariesLocal = 0; 1809 pcbddc->user_provided_isfordofs = PETSC_FALSE; 1810 pcbddc->n_ISForDofs = 0; 1811 pcbddc->n_ISForDofsLocal = 0; 1812 pcbddc->ISForDofs = 0; 1813 pcbddc->ISForDofsLocal = 0; 1814 pcbddc->ConstraintMatrix = 0; 1815 pcbddc->use_exact_dirichlet_trick = PETSC_TRUE; 1816 pcbddc->coarse_loc_to_glob = 0; 1817 pcbddc->coarsening_ratio = 8; 1818 pcbddc->current_level = 0; 1819 pcbddc->max_levels = 0; 1820 pcbddc->use_coarse_estimates = PETSC_FALSE; 1821 pcbddc->coarse_subassembling = 0; 1822 pcbddc->coarse_subassembling_init = 0; 1823 1824 /* create local graph structure */ 1825 ierr = PCBDDCGraphCreate(&pcbddc->mat_graph);CHKERRQ(ierr); 1826 1827 /* scaling */ 1828 pcbddc->work_scaling = 0; 1829 pcbddc->use_deluxe_scaling = PETSC_FALSE; 1830 pcbddc->deluxe_threshold = 1; 1831 pcbddc->deluxe_rebuild = PETSC_FALSE; 1832 pcbddc->deluxe_layers = -1; 1833 pcbddc->deluxe_use_useradj = PETSC_FALSE; 1834 pcbddc->deluxe_compute_rowadj = PETSC_TRUE; 1835 1836 /* function pointers */ 1837 pc->ops->apply = PCApply_BDDC; 1838 pc->ops->applytranspose = PCApplyTranspose_BDDC; 1839 pc->ops->setup = PCSetUp_BDDC; 1840 pc->ops->destroy = PCDestroy_BDDC; 1841 pc->ops->setfromoptions = PCSetFromOptions_BDDC; 1842 pc->ops->view = 0; 1843 pc->ops->applyrichardson = 0; 1844 pc->ops->applysymmetricleft = 0; 1845 pc->ops->applysymmetricright = 0; 1846 pc->ops->presolve = PCPreSolve_BDDC; 1847 pc->ops->postsolve = PCPostSolve_BDDC; 1848 1849 /* composing function */ 1850 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetChangeOfBasisLocalMat_C",PCBDDCSetChangeOfBasisLocalMat_BDDC);CHKERRQ(ierr); 1851 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetPrimalVerticesLocalIS_C",PCBDDCSetPrimalVerticesLocalIS_BDDC);CHKERRQ(ierr); 1852 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetCoarseningRatio_C",PCBDDCSetCoarseningRatio_BDDC);CHKERRQ(ierr); 1853 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetLevel_C",PCBDDCSetLevel_BDDC);CHKERRQ(ierr); 1854 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetUseExactDirichlet_C",PCBDDCSetUseExactDirichlet_BDDC);CHKERRQ(ierr); 1855 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetLevels_C",PCBDDCSetLevels_BDDC);CHKERRQ(ierr); 1856 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetNullSpace_C",PCBDDCSetNullSpace_BDDC);CHKERRQ(ierr); 1857 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetDirichletBoundaries_C",PCBDDCSetDirichletBoundaries_BDDC);CHKERRQ(ierr); 1858 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetDirichletBoundariesLocal_C",PCBDDCSetDirichletBoundariesLocal_BDDC);CHKERRQ(ierr); 1859 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetNeumannBoundaries_C",PCBDDCSetNeumannBoundaries_BDDC);CHKERRQ(ierr); 1860 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetNeumannBoundariesLocal_C",PCBDDCSetNeumannBoundariesLocal_BDDC);CHKERRQ(ierr); 1861 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCGetDirichletBoundaries_C",PCBDDCGetDirichletBoundaries_BDDC);CHKERRQ(ierr); 1862 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCGetDirichletBoundariesLocal_C",PCBDDCGetDirichletBoundariesLocal_BDDC);CHKERRQ(ierr); 1863 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCGetNeumannBoundaries_C",PCBDDCGetNeumannBoundaries_BDDC);CHKERRQ(ierr); 1864 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCGetNeumannBoundariesLocal_C",PCBDDCGetNeumannBoundariesLocal_BDDC);CHKERRQ(ierr); 1865 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetDofsSplitting_C",PCBDDCSetDofsSplitting_BDDC);CHKERRQ(ierr); 1866 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetDofsSplittingLocal_C",PCBDDCSetDofsSplittingLocal_BDDC);CHKERRQ(ierr); 1867 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetLocalAdjacencyGraph_C",PCBDDCSetLocalAdjacencyGraph_BDDC);CHKERRQ(ierr); 1868 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCCreateFETIDPOperators_C",PCBDDCCreateFETIDPOperators_BDDC);CHKERRQ(ierr); 1869 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCMatFETIDPGetRHS_C",PCBDDCMatFETIDPGetRHS_BDDC);CHKERRQ(ierr); 1870 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCMatFETIDPGetSolution_C",PCBDDCMatFETIDPGetSolution_BDDC);CHKERRQ(ierr); 1871 PetscFunctionReturn(0); 1872 } 1873 1874