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