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 ierr = PetscMalloc(n_is*sizeof(IS),&pcbddc->ISForDofsLocal);CHKERRQ(ierr); 804 for (i=0;i<n_is;i++) { 805 ierr = PetscObjectReference((PetscObject)ISForDofs[i]);CHKERRQ(ierr); 806 pcbddc->ISForDofsLocal[i]=ISForDofs[i]; 807 } 808 pcbddc->n_ISForDofsLocal=n_is; 809 if (n_is) pcbddc->user_provided_isfordofs = PETSC_TRUE; 810 pcbddc->recompute_topography = PETSC_TRUE; 811 PetscFunctionReturn(0); 812 } 813 814 #undef __FUNCT__ 815 #define __FUNCT__ "PCBDDCSetDofsSplittingLocal" 816 /*@ 817 PCBDDCSetDofsSplittingLocal - Set index sets defining fields of the local subdomain matrix 818 819 Collective 820 821 Input Parameters: 822 + pc - the preconditioning context 823 - n_is - number of index sets defining the fields 824 . ISForDofs - array of IS describing the fields in local ordering 825 826 Level: intermediate 827 828 Notes: n_is should be the same among processes. Not all nodes need to be listed: unlisted nodes will belong to a different field. 829 830 .seealso: PCBDDC 831 @*/ 832 PetscErrorCode PCBDDCSetDofsSplittingLocal(PC pc,PetscInt n_is, IS ISForDofs[]) 833 { 834 PetscInt i; 835 PetscErrorCode ierr; 836 837 PetscFunctionBegin; 838 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 839 PetscValidLogicalCollectiveInt(pc,n_is,2); 840 for (i=0;i<n_is;i++) { 841 PetscCheckSameComm(pc,1,ISForDofs[i],3); 842 PetscValidHeaderSpecific(ISForDofs[i],IS_CLASSID,3); 843 } 844 ierr = PetscTryMethod(pc,"PCBDDCSetDofsSplittingLocal_C",(PC,PetscInt,IS[]),(pc,n_is,ISForDofs));CHKERRQ(ierr); 845 PetscFunctionReturn(0); 846 } 847 /* -------------------------------------------------------------------------- */ 848 849 #undef __FUNCT__ 850 #define __FUNCT__ "PCBDDCSetDofsSplitting_BDDC" 851 static PetscErrorCode PCBDDCSetDofsSplitting_BDDC(PC pc,PetscInt n_is, IS ISForDofs[]) 852 { 853 PC_BDDC *pcbddc = (PC_BDDC*)pc->data; 854 PetscInt i; 855 PetscErrorCode ierr; 856 857 PetscFunctionBegin; 858 /* Destroy ISes if they were already set */ 859 for (i=0;i<pcbddc->n_ISForDofs;i++) { 860 ierr = ISDestroy(&pcbddc->ISForDofs[i]);CHKERRQ(ierr); 861 } 862 ierr = PetscFree(pcbddc->ISForDofs);CHKERRQ(ierr); 863 /* last user setting takes precendence -> destroy any other customization */ 864 for (i=0;i<pcbddc->n_ISForDofsLocal;i++) { 865 ierr = ISDestroy(&pcbddc->ISForDofsLocal[i]);CHKERRQ(ierr); 866 } 867 ierr = PetscFree(pcbddc->ISForDofsLocal);CHKERRQ(ierr); 868 pcbddc->n_ISForDofsLocal = 0; 869 /* allocate space then set */ 870 ierr = PetscMalloc1(n_is,&pcbddc->ISForDofs);CHKERRQ(ierr); 871 for (i=0;i<n_is;i++) { 872 ierr = PetscObjectReference((PetscObject)ISForDofs[i]);CHKERRQ(ierr); 873 pcbddc->ISForDofs[i]=ISForDofs[i]; 874 } 875 pcbddc->n_ISForDofs=n_is; 876 if (n_is) pcbddc->user_provided_isfordofs = PETSC_TRUE; 877 pcbddc->recompute_topography = PETSC_TRUE; 878 PetscFunctionReturn(0); 879 } 880 881 #undef __FUNCT__ 882 #define __FUNCT__ "PCBDDCSetDofsSplitting" 883 /*@ 884 PCBDDCSetDofsSplitting - Set index sets defining fields of the global matrix 885 886 Collective 887 888 Input Parameters: 889 + pc - the preconditioning context 890 - n_is - number of index sets defining the fields 891 . ISForDofs - array of IS describing the fields in global ordering 892 893 Level: intermediate 894 895 Notes: Any process can list any global node. Not all nodes need to be listed: unlisted nodes will belong to a different field. 896 897 .seealso: PCBDDC 898 @*/ 899 PetscErrorCode PCBDDCSetDofsSplitting(PC pc,PetscInt n_is, IS ISForDofs[]) 900 { 901 PetscInt i; 902 PetscErrorCode ierr; 903 904 PetscFunctionBegin; 905 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 906 PetscValidLogicalCollectiveInt(pc,n_is,2); 907 for (i=0;i<n_is;i++) { 908 PetscCheckSameComm(pc,1,ISForDofs[i],3); 909 PetscValidHeaderSpecific(ISForDofs[i],IS_CLASSID,3); 910 } 911 ierr = PetscTryMethod(pc,"PCBDDCSetDofsSplitting_C",(PC,PetscInt,IS[]),(pc,n_is,ISForDofs));CHKERRQ(ierr); 912 PetscFunctionReturn(0); 913 } 914 915 /* -------------------------------------------------------------------------- */ 916 #undef __FUNCT__ 917 #define __FUNCT__ "PCPreSolve_BDDC" 918 /* -------------------------------------------------------------------------- */ 919 /* 920 PCPreSolve_BDDC - Changes the right hand side and (if necessary) the initial 921 guess if a transformation of basis approach has been selected. 922 923 Input Parameter: 924 + pc - the preconditioner contex 925 926 Application Interface Routine: PCPreSolve() 927 928 Notes: 929 The interface routine PCPreSolve() is not usually called directly by 930 the user, but instead is called by KSPSolve(). 931 */ 932 static PetscErrorCode PCPreSolve_BDDC(PC pc, KSP ksp, Vec rhs, Vec x) 933 { 934 PetscErrorCode ierr; 935 PC_BDDC *pcbddc = (PC_BDDC*)pc->data; 936 PC_IS *pcis = (PC_IS*)(pc->data); 937 IS dirIS; 938 Vec used_vec; 939 PetscBool copy_rhs = PETSC_TRUE; 940 941 PetscFunctionBegin; 942 /* if we are working with cg, one dirichlet solve can be avoided during Krylov iterations */ 943 if (ksp) { 944 PetscBool iscg; 945 ierr = PetscObjectTypeCompare((PetscObject)ksp,KSPCG,&iscg);CHKERRQ(ierr); 946 if (!iscg) { 947 ierr = PCBDDCSetUseExactDirichlet(pc,PETSC_FALSE);CHKERRQ(ierr); 948 } 949 } 950 /* Creates parallel work vectors used in presolve */ 951 if (!pcbddc->original_rhs) { 952 ierr = VecDuplicate(pcis->vec1_global,&pcbddc->original_rhs);CHKERRQ(ierr); 953 } 954 if (!pcbddc->temp_solution) { 955 ierr = VecDuplicate(pcis->vec1_global,&pcbddc->temp_solution);CHKERRQ(ierr); 956 } 957 958 if (x) { 959 ierr = PetscObjectReference((PetscObject)x);CHKERRQ(ierr); 960 used_vec = x; 961 } else { /* it can only happen when calling PCBDDCMatFETIDPGetRHS */ 962 ierr = PetscObjectReference((PetscObject)pcbddc->temp_solution);CHKERRQ(ierr); 963 used_vec = pcbddc->temp_solution; 964 ierr = VecSet(used_vec,0.0);CHKERRQ(ierr); 965 } 966 967 /* hack into ksp data structure since PCPreSolve comes earlier than setting to zero the guess in src/ksp/ksp/interface/itfunc.c */ 968 if (ksp) { 969 ierr = KSPGetInitialGuessNonzero(ksp,&pcbddc->ksp_guess_nonzero);CHKERRQ(ierr); 970 if (!pcbddc->ksp_guess_nonzero) { 971 ierr = VecSet(used_vec,0.0);CHKERRQ(ierr); 972 } 973 } 974 975 pcbddc->rhs_change = PETSC_FALSE; 976 977 /* Take into account zeroed rows -> change rhs and store solution removed */ 978 /* note that Dirichlet boundaries in global ordering (if any) has already been translated into local ordering in PCBDDCAnalyzeInterface */ 979 ierr = PCBDDCGetDirichletBoundariesLocal(pc,&dirIS);CHKERRQ(ierr); 980 if (rhs && dirIS) { 981 Mat_IS *matis = (Mat_IS*)pc->pmat->data; 982 PetscInt dirsize,i,*is_indices; 983 PetscScalar *array_x,*array_diagonal; 984 985 ierr = MatGetDiagonal(pc->pmat,pcis->vec1_global);CHKERRQ(ierr); 986 ierr = VecPointwiseDivide(pcis->vec1_global,rhs,pcis->vec1_global);CHKERRQ(ierr); 987 ierr = VecScatterBegin(matis->ctx,pcis->vec1_global,pcis->vec2_N,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 988 ierr = VecScatterEnd(matis->ctx,pcis->vec1_global,pcis->vec2_N,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 989 ierr = VecScatterBegin(matis->ctx,used_vec,pcis->vec1_N,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 990 ierr = VecScatterEnd(matis->ctx,used_vec,pcis->vec1_N,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 991 ierr = ISGetLocalSize(dirIS,&dirsize);CHKERRQ(ierr); 992 ierr = VecGetArray(pcis->vec1_N,&array_x);CHKERRQ(ierr); 993 ierr = VecGetArray(pcis->vec2_N,&array_diagonal);CHKERRQ(ierr); 994 ierr = ISGetIndices(dirIS,(const PetscInt**)&is_indices);CHKERRQ(ierr); 995 for (i=0; i<dirsize; i++) array_x[is_indices[i]] = array_diagonal[is_indices[i]]; 996 ierr = ISRestoreIndices(dirIS,(const PetscInt**)&is_indices);CHKERRQ(ierr); 997 ierr = VecRestoreArray(pcis->vec2_N,&array_diagonal);CHKERRQ(ierr); 998 ierr = VecRestoreArray(pcis->vec1_N,&array_x);CHKERRQ(ierr); 999 ierr = VecScatterBegin(matis->ctx,pcis->vec1_N,used_vec,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1000 ierr = VecScatterEnd(matis->ctx,pcis->vec1_N,used_vec,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1001 pcbddc->rhs_change = PETSC_TRUE; 1002 } 1003 1004 /* remove the computed solution or the initial guess from the rhs */ 1005 if (pcbddc->rhs_change || (ksp && pcbddc->ksp_guess_nonzero) ) { 1006 /* store the original rhs */ 1007 if (copy_rhs) { 1008 ierr = VecCopy(rhs,pcbddc->original_rhs);CHKERRQ(ierr); 1009 copy_rhs = PETSC_FALSE; 1010 } 1011 pcbddc->rhs_change = PETSC_TRUE; 1012 ierr = VecScale(used_vec,-1.0);CHKERRQ(ierr); 1013 ierr = MatMultAdd(pc->pmat,used_vec,rhs,rhs);CHKERRQ(ierr); 1014 ierr = VecScale(used_vec,-1.0);CHKERRQ(ierr); 1015 ierr = VecCopy(used_vec,pcbddc->temp_solution);CHKERRQ(ierr); 1016 } 1017 ierr = VecDestroy(&used_vec);CHKERRQ(ierr); 1018 1019 /* store partially computed solution and set initial guess */ 1020 if (x) { 1021 ierr = VecSet(x,0.0);CHKERRQ(ierr); 1022 if (pcbddc->use_exact_dirichlet_trick) { 1023 ierr = VecScatterBegin(pcis->global_to_D,rhs,pcis->vec1_D,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1024 ierr = VecScatterEnd(pcis->global_to_D,rhs,pcis->vec1_D,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1025 ierr = KSPSolve(pcbddc->ksp_D,pcis->vec1_D,pcis->vec2_D);CHKERRQ(ierr); 1026 ierr = VecScatterBegin(pcis->global_to_D,pcis->vec2_D,x,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1027 ierr = VecScatterEnd(pcis->global_to_D,pcis->vec2_D,x,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1028 if (ksp && !pcbddc->ksp_guess_nonzero) { 1029 ierr = KSPSetInitialGuessNonzero(ksp,PETSC_TRUE);CHKERRQ(ierr); 1030 } 1031 } 1032 } 1033 1034 if (pcbddc->ChangeOfBasisMatrix) { 1035 PCBDDCChange_ctx change_ctx; 1036 1037 /* get change ctx */ 1038 ierr = MatShellGetContext(pcbddc->new_global_mat,&change_ctx);CHKERRQ(ierr); 1039 1040 /* set current iteration matrix inside change context (change of basis has been already set into the ctx during PCSetUp) */ 1041 ierr = MatDestroy(&change_ctx->original_mat);CHKERRQ(ierr); 1042 ierr = PetscObjectReference((PetscObject)pc->mat);CHKERRQ(ierr); 1043 change_ctx->original_mat = pc->mat; 1044 1045 /* change iteration matrix */ 1046 ierr = MatDestroy(&pc->mat);CHKERRQ(ierr); 1047 ierr = PetscObjectReference((PetscObject)pcbddc->new_global_mat);CHKERRQ(ierr); 1048 pc->mat = pcbddc->new_global_mat; 1049 1050 /* store the original rhs */ 1051 if (copy_rhs) { 1052 ierr = VecCopy(rhs,pcbddc->original_rhs);CHKERRQ(ierr); 1053 copy_rhs = PETSC_FALSE; 1054 } 1055 1056 /* change rhs */ 1057 ierr = MatMultTranspose(change_ctx->global_change,rhs,pcis->vec1_global);CHKERRQ(ierr); 1058 ierr = VecCopy(pcis->vec1_global,rhs);CHKERRQ(ierr); 1059 pcbddc->rhs_change = PETSC_TRUE; 1060 } 1061 1062 /* remove nullspace if present */ 1063 if (ksp && x && pcbddc->NullSpace) { 1064 ierr = MatNullSpaceRemove(pcbddc->NullSpace,x);CHKERRQ(ierr); 1065 /* store the original rhs */ 1066 if (copy_rhs) { 1067 ierr = VecCopy(rhs,pcbddc->original_rhs);CHKERRQ(ierr); 1068 copy_rhs = PETSC_FALSE; 1069 } 1070 pcbddc->rhs_change = PETSC_TRUE; 1071 ierr = MatNullSpaceRemove(pcbddc->NullSpace,rhs);CHKERRQ(ierr); 1072 } 1073 PetscFunctionReturn(0); 1074 } 1075 1076 /* -------------------------------------------------------------------------- */ 1077 #undef __FUNCT__ 1078 #define __FUNCT__ "PCPostSolve_BDDC" 1079 /* -------------------------------------------------------------------------- */ 1080 /* 1081 PCPostSolve_BDDC - Changes the computed solution if a transformation of basis 1082 approach has been selected. Also, restores rhs to its original state. 1083 1084 Input Parameter: 1085 + pc - the preconditioner contex 1086 1087 Application Interface Routine: PCPostSolve() 1088 1089 Notes: 1090 The interface routine PCPostSolve() is not usually called directly by 1091 the user, but instead is called by KSPSolve(). 1092 */ 1093 static PetscErrorCode PCPostSolve_BDDC(PC pc, KSP ksp, Vec rhs, Vec x) 1094 { 1095 PetscErrorCode ierr; 1096 PC_BDDC *pcbddc = (PC_BDDC*)pc->data; 1097 1098 PetscFunctionBegin; 1099 if (pcbddc->ChangeOfBasisMatrix) { 1100 PCBDDCChange_ctx change_ctx; 1101 1102 /* get change ctx */ 1103 ierr = MatShellGetContext(pcbddc->new_global_mat,&change_ctx);CHKERRQ(ierr); 1104 1105 /* restore iteration matrix */ 1106 ierr = MatDestroy(&pc->mat);CHKERRQ(ierr); 1107 ierr = PetscObjectReference((PetscObject)change_ctx->original_mat);CHKERRQ(ierr); 1108 pc->mat = change_ctx->original_mat; 1109 1110 /* get solution in original basis */ 1111 if (x) { 1112 PC_IS *pcis = (PC_IS*)(pc->data); 1113 ierr = MatMult(change_ctx->global_change,x,pcis->vec1_global);CHKERRQ(ierr); 1114 ierr = VecCopy(pcis->vec1_global,x);CHKERRQ(ierr); 1115 } 1116 } 1117 1118 /* add solution removed in presolve */ 1119 if (x && pcbddc->rhs_change) { 1120 ierr = VecAXPY(x,1.0,pcbddc->temp_solution);CHKERRQ(ierr); 1121 } 1122 1123 /* restore rhs to its original state */ 1124 if (rhs && pcbddc->rhs_change) { 1125 ierr = VecCopy(pcbddc->original_rhs,rhs);CHKERRQ(ierr); 1126 } 1127 pcbddc->rhs_change = PETSC_FALSE; 1128 1129 /* restore ksp guess state */ 1130 if (ksp) { 1131 ierr = KSPSetInitialGuessNonzero(ksp,pcbddc->ksp_guess_nonzero);CHKERRQ(ierr); 1132 } 1133 PetscFunctionReturn(0); 1134 } 1135 /* -------------------------------------------------------------------------- */ 1136 #undef __FUNCT__ 1137 #define __FUNCT__ "PCSetUp_BDDC" 1138 /* -------------------------------------------------------------------------- */ 1139 /* 1140 PCSetUp_BDDC - Prepares for the use of the BDDC preconditioner 1141 by setting data structures and options. 1142 1143 Input Parameter: 1144 + pc - the preconditioner context 1145 1146 Application Interface Routine: PCSetUp() 1147 1148 Notes: 1149 The interface routine PCSetUp() is not usually called directly by 1150 the user, but instead is called by PCApply() if necessary. 1151 */ 1152 PetscErrorCode PCSetUp_BDDC(PC pc) 1153 { 1154 PetscErrorCode ierr; 1155 PC_BDDC* pcbddc = (PC_BDDC*)pc->data; 1156 MatNullSpace nearnullspace; 1157 PetscBool computeis,computetopography,computesolvers; 1158 PetscBool new_nearnullspace_provided; 1159 1160 PetscFunctionBegin; 1161 /* the following lines of code should be replaced by a better logic between PCIS, PCNN, PCBDDC and other future nonoverlapping preconditioners */ 1162 /* For BDDC we need to define a local "Neumann" problem different to that defined in PCISSetup 1163 Also, BDDC directly build the Dirichlet problem */ 1164 /* split work */ 1165 if (pc->setupcalled) { 1166 computeis = PETSC_FALSE; 1167 if (pc->flag == SAME_NONZERO_PATTERN) { 1168 computetopography = PETSC_FALSE; 1169 computesolvers = PETSC_TRUE; 1170 } else { /* DIFFERENT_NONZERO_PATTERN */ 1171 computetopography = PETSC_TRUE; 1172 computesolvers = PETSC_TRUE; 1173 } 1174 } else { 1175 computeis = PETSC_TRUE; 1176 computetopography = PETSC_TRUE; 1177 computesolvers = PETSC_TRUE; 1178 } 1179 if (pcbddc->recompute_topography) { 1180 computetopography = PETSC_TRUE; 1181 } 1182 1183 /* Get stdout for dbg */ 1184 if (pcbddc->dbg_flag) { 1185 if (!pcbddc->dbg_viewer) { 1186 pcbddc->dbg_viewer = PETSC_VIEWER_STDOUT_(PetscObjectComm((PetscObject)pc)); 1187 ierr = PetscViewerASCIISynchronizedAllow(pcbddc->dbg_viewer,PETSC_TRUE);CHKERRQ(ierr); 1188 } 1189 ierr = PetscViewerASCIIAddTab(pcbddc->dbg_viewer,2*pcbddc->current_level);CHKERRQ(ierr); 1190 } 1191 1192 /* Set up all the "iterative substructuring" common block without computing solvers */ 1193 if (computeis) { 1194 /* HACK INTO PCIS */ 1195 PC_IS* pcis = (PC_IS*)pc->data; 1196 pcis->computesolvers = PETSC_FALSE; 1197 ierr = PCISSetUp(pc);CHKERRQ(ierr); 1198 ierr = ISLocalToGlobalMappingCreateIS(pcis->is_B_local,&pcbddc->BtoNmap);CHKERRQ(ierr); 1199 } 1200 1201 /* use_change_of_basis flag is used to automatically compute a change of basis from constraints */ 1202 if (pcbddc->user_ChangeOfBasisMatrix) { 1203 pcbddc->use_change_of_basis = PETSC_FALSE; 1204 } 1205 1206 /* Analyze interface */ 1207 if (computetopography) { 1208 ierr = PCBDDCAnalyzeInterface(pc);CHKERRQ(ierr); 1209 /* Schurs on subsets should be reset */ 1210 if (pcbddc->deluxe_ctx) { 1211 ierr = PCBDDCSubSchursReset(pcbddc->deluxe_ctx->sub_schurs);CHKERRQ(ierr); 1212 } 1213 } 1214 1215 /* infer if NullSpace object attached to Mat via MatSetNearNullSpace has changed */ 1216 new_nearnullspace_provided = PETSC_FALSE; 1217 ierr = MatGetNearNullSpace(pc->pmat,&nearnullspace);CHKERRQ(ierr); 1218 if (pcbddc->onearnullspace) { /* already used nearnullspace */ 1219 if (!nearnullspace) { /* near null space attached to mat has been destroyed */ 1220 new_nearnullspace_provided = PETSC_TRUE; 1221 } else { 1222 /* determine if the two nullspaces are different (should be lightweight) */ 1223 if (nearnullspace != pcbddc->onearnullspace) { 1224 new_nearnullspace_provided = PETSC_TRUE; 1225 } else { /* maybe the user has changed the content of the nearnullspace so check vectors ObjectStateId */ 1226 PetscInt i; 1227 const Vec *nearnullvecs; 1228 PetscObjectState state; 1229 PetscInt nnsp_size; 1230 ierr = MatNullSpaceGetVecs(nearnullspace,NULL,&nnsp_size,&nearnullvecs);CHKERRQ(ierr); 1231 for (i=0;i<nnsp_size;i++) { 1232 ierr = PetscObjectStateGet((PetscObject)nearnullvecs[i],&state);CHKERRQ(ierr); 1233 if (pcbddc->onearnullvecs_state[i] != state) { 1234 new_nearnullspace_provided = PETSC_TRUE; 1235 break; 1236 } 1237 } 1238 } 1239 } 1240 } else { 1241 if (!nearnullspace) { /* both nearnullspaces are null */ 1242 new_nearnullspace_provided = PETSC_FALSE; 1243 } else { /* nearnullspace attached later */ 1244 new_nearnullspace_provided = PETSC_TRUE; 1245 } 1246 } 1247 1248 /* Setup constraints and related work vectors */ 1249 /* reset primal space flags */ 1250 pcbddc->new_primal_space = PETSC_FALSE; 1251 pcbddc->new_primal_space_local = PETSC_FALSE; 1252 if (computetopography || new_nearnullspace_provided) { 1253 /* It also sets the primal space flags */ 1254 ierr = PCBDDCConstraintsSetUp(pc);CHKERRQ(ierr); 1255 /* Allocate needed local vectors (which depends on quantities defined during ConstraintsSetUp) */ 1256 ierr = PCBDDCSetUpLocalWorkVectors(pc);CHKERRQ(ierr); 1257 /* Schurs on subsets should be reset */ 1258 if (pcbddc->deluxe_ctx) { 1259 ierr = PCBDDCSubSchursReset(pcbddc->deluxe_ctx->sub_schurs);CHKERRQ(ierr); 1260 } 1261 } 1262 1263 if (computesolvers || pcbddc->new_primal_space) { 1264 /* Create coarse and local stuffs */ 1265 ierr = PCBDDCSetUpSolvers(pc);CHKERRQ(ierr); 1266 /* Create scaling operators */ 1267 ierr = PCBDDCScalingSetUp(pc);CHKERRQ(ierr); 1268 } 1269 1270 if (pcbddc->dbg_flag) { 1271 ierr = PetscViewerASCIISubtractTab(pcbddc->dbg_viewer,2*pcbddc->current_level);CHKERRQ(ierr); 1272 } 1273 PetscFunctionReturn(0); 1274 } 1275 1276 /* -------------------------------------------------------------------------- */ 1277 /* 1278 PCApply_BDDC - Applies the BDDC operator to a vector. 1279 1280 Input Parameters: 1281 . pc - the preconditioner context 1282 . r - input vector (global) 1283 1284 Output Parameter: 1285 . z - output vector (global) 1286 1287 Application Interface Routine: PCApply() 1288 */ 1289 #undef __FUNCT__ 1290 #define __FUNCT__ "PCApply_BDDC" 1291 PetscErrorCode PCApply_BDDC(PC pc,Vec r,Vec z) 1292 { 1293 PC_IS *pcis = (PC_IS*)(pc->data); 1294 PC_BDDC *pcbddc = (PC_BDDC*)(pc->data); 1295 PetscErrorCode ierr; 1296 const PetscScalar one = 1.0; 1297 const PetscScalar m_one = -1.0; 1298 const PetscScalar zero = 0.0; 1299 1300 /* This code is similar to that provided in nn.c for PCNN 1301 NN interface preconditioner changed to BDDC 1302 Added support for M_3 preconditioner in the reference article (code is active if pcbddc->switch_static = PETSC_TRUE) */ 1303 1304 PetscFunctionBegin; 1305 if (!pcbddc->use_exact_dirichlet_trick) { 1306 /* First Dirichlet solve */ 1307 ierr = VecScatterBegin(pcis->global_to_D,r,pcis->vec1_D,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1308 ierr = VecScatterEnd(pcis->global_to_D,r,pcis->vec1_D,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1309 ierr = KSPSolve(pcbddc->ksp_D,pcis->vec1_D,pcis->vec2_D);CHKERRQ(ierr); 1310 /* 1311 Assembling right hand side for BDDC operator 1312 - pcis->vec1_D for the Dirichlet part (if needed, i.e. prec_flag=PETSC_TRUE) 1313 - pcis->vec1_B the interface part of the global vector z 1314 */ 1315 ierr = VecScale(pcis->vec2_D,m_one);CHKERRQ(ierr); 1316 ierr = MatMult(pcis->A_BI,pcis->vec2_D,pcis->vec1_B);CHKERRQ(ierr); 1317 if (pcbddc->switch_static) { ierr = MatMultAdd(pcis->A_II,pcis->vec2_D,pcis->vec1_D,pcis->vec1_D);CHKERRQ(ierr); } 1318 ierr = VecScale(pcis->vec2_D,m_one);CHKERRQ(ierr); 1319 ierr = VecCopy(r,z);CHKERRQ(ierr); 1320 ierr = VecScatterBegin(pcis->global_to_B,pcis->vec1_B,z,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1321 ierr = VecScatterEnd(pcis->global_to_B,pcis->vec1_B,z,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1322 ierr = PCBDDCScalingRestriction(pc,z,pcis->vec1_B);CHKERRQ(ierr); 1323 } else { 1324 ierr = VecSet(pcis->vec1_D,zero);CHKERRQ(ierr); 1325 ierr = VecSet(pcis->vec2_D,zero);CHKERRQ(ierr); 1326 ierr = PCBDDCScalingRestriction(pc,r,pcis->vec1_B);CHKERRQ(ierr); 1327 } 1328 1329 /* Apply interface preconditioner 1330 input/output vecs: pcis->vec1_B and pcis->vec1_D */ 1331 ierr = PCBDDCApplyInterfacePreconditioner(pc,PETSC_FALSE);CHKERRQ(ierr); 1332 1333 /* Apply transpose of partition of unity operator */ 1334 ierr = PCBDDCScalingExtension(pc,pcis->vec1_B,z);CHKERRQ(ierr); 1335 1336 /* Second Dirichlet solve and assembling of output */ 1337 ierr = VecScatterBegin(pcis->global_to_B,z,pcis->vec1_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1338 ierr = VecScatterEnd(pcis->global_to_B,z,pcis->vec1_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1339 ierr = MatMult(pcis->A_IB,pcis->vec1_B,pcis->vec3_D);CHKERRQ(ierr); 1340 if (pcbddc->switch_static) { ierr = MatMultAdd(pcis->A_II,pcis->vec1_D,pcis->vec3_D,pcis->vec3_D);CHKERRQ(ierr); } 1341 ierr = KSPSolve(pcbddc->ksp_D,pcis->vec3_D,pcis->vec4_D);CHKERRQ(ierr); 1342 ierr = VecScale(pcis->vec4_D,m_one);CHKERRQ(ierr); 1343 if (pcbddc->switch_static) { ierr = VecAXPY(pcis->vec4_D,one,pcis->vec1_D);CHKERRQ(ierr); } 1344 ierr = VecAXPY(pcis->vec2_D,one,pcis->vec4_D);CHKERRQ(ierr); 1345 ierr = VecScatterBegin(pcis->global_to_D,pcis->vec2_D,z,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1346 ierr = VecScatterEnd(pcis->global_to_D,pcis->vec2_D,z,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1347 PetscFunctionReturn(0); 1348 } 1349 1350 /* -------------------------------------------------------------------------- */ 1351 /* 1352 PCApplyTranspose_BDDC - Applies the transpose of the BDDC operator to a vector. 1353 1354 Input Parameters: 1355 . pc - the preconditioner context 1356 . r - input vector (global) 1357 1358 Output Parameter: 1359 . z - output vector (global) 1360 1361 Application Interface Routine: PCApplyTranspose() 1362 */ 1363 #undef __FUNCT__ 1364 #define __FUNCT__ "PCApplyTranspose_BDDC" 1365 PetscErrorCode PCApplyTranspose_BDDC(PC pc,Vec r,Vec z) 1366 { 1367 PC_IS *pcis = (PC_IS*)(pc->data); 1368 PC_BDDC *pcbddc = (PC_BDDC*)(pc->data); 1369 PetscErrorCode ierr; 1370 const PetscScalar one = 1.0; 1371 const PetscScalar m_one = -1.0; 1372 const PetscScalar zero = 0.0; 1373 1374 PetscFunctionBegin; 1375 if (!pcbddc->use_exact_dirichlet_trick) { 1376 /* First Dirichlet solve */ 1377 ierr = VecScatterBegin(pcis->global_to_D,r,pcis->vec1_D,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1378 ierr = VecScatterEnd(pcis->global_to_D,r,pcis->vec1_D,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1379 ierr = KSPSolveTranspose(pcbddc->ksp_D,pcis->vec1_D,pcis->vec2_D);CHKERRQ(ierr); 1380 /* 1381 Assembling right hand side for BDDC operator 1382 - pcis->vec1_D for the Dirichlet part (if needed, i.e. prec_flag=PETSC_TRUE) 1383 - pcis->vec1_B the interface part of the global vector z 1384 */ 1385 ierr = VecScale(pcis->vec2_D,m_one);CHKERRQ(ierr); 1386 ierr = MatMultTranspose(pcis->A_IB,pcis->vec2_D,pcis->vec1_B);CHKERRQ(ierr); 1387 if (pcbddc->switch_static) { ierr = MatMultTransposeAdd(pcis->A_II,pcis->vec2_D,pcis->vec1_D,pcis->vec1_D);CHKERRQ(ierr); } 1388 ierr = VecScale(pcis->vec2_D,m_one);CHKERRQ(ierr); 1389 ierr = VecCopy(r,z);CHKERRQ(ierr); 1390 ierr = VecScatterBegin(pcis->global_to_B,pcis->vec1_B,z,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1391 ierr = VecScatterEnd(pcis->global_to_B,pcis->vec1_B,z,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1392 ierr = PCBDDCScalingRestriction(pc,z,pcis->vec1_B);CHKERRQ(ierr); 1393 } else { 1394 ierr = VecSet(pcis->vec1_D,zero);CHKERRQ(ierr); 1395 ierr = VecSet(pcis->vec2_D,zero);CHKERRQ(ierr); 1396 ierr = PCBDDCScalingRestriction(pc,r,pcis->vec1_B);CHKERRQ(ierr); 1397 } 1398 1399 /* Apply interface preconditioner 1400 input/output vecs: pcis->vec1_B and pcis->vec1_D */ 1401 ierr = PCBDDCApplyInterfacePreconditioner(pc,PETSC_TRUE);CHKERRQ(ierr); 1402 1403 /* Apply transpose of partition of unity operator */ 1404 ierr = PCBDDCScalingExtension(pc,pcis->vec1_B,z);CHKERRQ(ierr); 1405 1406 /* Second Dirichlet solve and assembling of output */ 1407 ierr = VecScatterBegin(pcis->global_to_B,z,pcis->vec1_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1408 ierr = VecScatterEnd(pcis->global_to_B,z,pcis->vec1_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1409 ierr = MatMultTranspose(pcis->A_BI,pcis->vec1_B,pcis->vec3_D);CHKERRQ(ierr); 1410 if (pcbddc->switch_static) { ierr = MatMultTransposeAdd(pcis->A_II,pcis->vec1_D,pcis->vec3_D,pcis->vec3_D);CHKERRQ(ierr); } 1411 ierr = KSPSolveTranspose(pcbddc->ksp_D,pcis->vec3_D,pcis->vec4_D);CHKERRQ(ierr); 1412 ierr = VecScale(pcis->vec4_D,m_one);CHKERRQ(ierr); 1413 if (pcbddc->switch_static) { ierr = VecAXPY(pcis->vec4_D,one,pcis->vec1_D);CHKERRQ(ierr); } 1414 ierr = VecAXPY(pcis->vec2_D,one,pcis->vec4_D);CHKERRQ(ierr); 1415 ierr = VecScatterBegin(pcis->global_to_D,pcis->vec2_D,z,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1416 ierr = VecScatterEnd(pcis->global_to_D,pcis->vec2_D,z,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1417 PetscFunctionReturn(0); 1418 } 1419 /* -------------------------------------------------------------------------- */ 1420 1421 #undef __FUNCT__ 1422 #define __FUNCT__ "PCDestroy_BDDC" 1423 PetscErrorCode PCDestroy_BDDC(PC pc) 1424 { 1425 PC_BDDC *pcbddc = (PC_BDDC*)pc->data; 1426 PetscErrorCode ierr; 1427 1428 PetscFunctionBegin; 1429 /* free data created by PCIS */ 1430 ierr = PCISDestroy(pc);CHKERRQ(ierr); 1431 /* free BDDC custom data */ 1432 ierr = PCBDDCResetCustomization(pc);CHKERRQ(ierr); 1433 /* destroy objects related to topography */ 1434 ierr = PCBDDCResetTopography(pc);CHKERRQ(ierr); 1435 /* free allocated graph structure */ 1436 ierr = PetscFree(pcbddc->mat_graph);CHKERRQ(ierr); 1437 /* destroy objects for scaling operator */ 1438 ierr = PCBDDCScalingDestroy(pc);CHKERRQ(ierr); 1439 ierr = PetscFree(pcbddc->deluxe_ctx);CHKERRQ(ierr); 1440 /* free solvers stuff */ 1441 ierr = PCBDDCResetSolvers(pc);CHKERRQ(ierr); 1442 /* free global vectors needed in presolve */ 1443 ierr = VecDestroy(&pcbddc->temp_solution);CHKERRQ(ierr); 1444 ierr = VecDestroy(&pcbddc->original_rhs);CHKERRQ(ierr); 1445 /* free stuff for change of basis hooks */ 1446 if (pcbddc->new_global_mat) { 1447 PCBDDCChange_ctx change_ctx; 1448 ierr = MatShellGetContext(pcbddc->new_global_mat,&change_ctx);CHKERRQ(ierr); 1449 ierr = MatDestroy(&change_ctx->original_mat);CHKERRQ(ierr); 1450 ierr = MatDestroy(&change_ctx->global_change);CHKERRQ(ierr); 1451 ierr = VecDestroyVecs(2,&change_ctx->work);CHKERRQ(ierr); 1452 ierr = PetscFree(change_ctx);CHKERRQ(ierr); 1453 } 1454 ierr = MatDestroy(&pcbddc->new_global_mat);CHKERRQ(ierr); 1455 /* remove map from local boundary to local numbering */ 1456 ierr = ISLocalToGlobalMappingDestroy(&pcbddc->BtoNmap);CHKERRQ(ierr); 1457 /* remove functions */ 1458 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetChangeOfBasisMat_C",NULL);CHKERRQ(ierr); 1459 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetPrimalVerticesLocalIS_C",NULL);CHKERRQ(ierr); 1460 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetCoarseningRatio_C",NULL);CHKERRQ(ierr); 1461 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetLevel_C",NULL);CHKERRQ(ierr); 1462 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetUseExactDirichlet_C",NULL);CHKERRQ(ierr); 1463 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetLevels_C",NULL);CHKERRQ(ierr); 1464 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetNullSpace_C",NULL);CHKERRQ(ierr); 1465 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetDirichletBoundaries_C",NULL);CHKERRQ(ierr); 1466 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetDirichletBoundariesLocal_C",NULL);CHKERRQ(ierr); 1467 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetNeumannBoundaries_C",NULL);CHKERRQ(ierr); 1468 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetNeumannBoundariesLocal_C",NULL);CHKERRQ(ierr); 1469 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCGetDirichletBoundaries_C",NULL);CHKERRQ(ierr); 1470 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCGetDirichletBoundariesLocal_C",NULL);CHKERRQ(ierr); 1471 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCGetNeumannBoundaries_C",NULL);CHKERRQ(ierr); 1472 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCGetNeumannBoundariesLocal_C",NULL);CHKERRQ(ierr); 1473 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetDofsSplitting_C",NULL);CHKERRQ(ierr); 1474 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetDofsSplittingLocal_C",NULL);CHKERRQ(ierr); 1475 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetLocalAdjacencyGraph_C",NULL);CHKERRQ(ierr); 1476 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCCreateFETIDPOperators_C",NULL);CHKERRQ(ierr); 1477 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCMatFETIDPGetRHS_C",NULL);CHKERRQ(ierr); 1478 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCMatFETIDPGetSolution_C",NULL);CHKERRQ(ierr); 1479 /* Free the private data structure */ 1480 ierr = PetscFree(pc->data);CHKERRQ(ierr); 1481 PetscFunctionReturn(0); 1482 } 1483 /* -------------------------------------------------------------------------- */ 1484 1485 #undef __FUNCT__ 1486 #define __FUNCT__ "PCBDDCMatFETIDPGetRHS_BDDC" 1487 static PetscErrorCode PCBDDCMatFETIDPGetRHS_BDDC(Mat fetidp_mat, Vec standard_rhs, Vec fetidp_flux_rhs) 1488 { 1489 FETIDPMat_ctx mat_ctx; 1490 PC_IS* pcis; 1491 PC_BDDC* pcbddc; 1492 PetscErrorCode ierr; 1493 1494 PetscFunctionBegin; 1495 ierr = MatShellGetContext(fetidp_mat,&mat_ctx);CHKERRQ(ierr); 1496 pcis = (PC_IS*)mat_ctx->pc->data; 1497 pcbddc = (PC_BDDC*)mat_ctx->pc->data; 1498 1499 /* change of basis for physical rhs if needed 1500 It also changes the rhs in case of dirichlet boundaries */ 1501 ierr = PCPreSolve_BDDC(mat_ctx->pc,NULL,standard_rhs,NULL);CHKERRQ(ierr); 1502 /* store vectors for computation of fetidp final solution */ 1503 ierr = VecScatterBegin(pcis->global_to_D,standard_rhs,mat_ctx->temp_solution_D,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1504 ierr = VecScatterEnd(pcis->global_to_D,standard_rhs,mat_ctx->temp_solution_D,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1505 /* scale rhs since it should be unassembled */ 1506 /* TODO use counter scaling? (also below) */ 1507 ierr = VecScatterBegin(pcis->global_to_B,standard_rhs,mat_ctx->temp_solution_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1508 ierr = VecScatterEnd(pcis->global_to_B,standard_rhs,mat_ctx->temp_solution_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1509 /* Apply partition of unity */ 1510 ierr = VecPointwiseMult(mat_ctx->temp_solution_B,pcis->D,mat_ctx->temp_solution_B);CHKERRQ(ierr); 1511 /* ierr = PCBDDCScalingRestriction(mat_ctx->pc,standard_rhs,mat_ctx->temp_solution_B);CHKERRQ(ierr); */ 1512 if (!pcbddc->switch_static) { 1513 /* compute partially subassembled Schur complement right-hand side */ 1514 ierr = KSPSolve(pcbddc->ksp_D,mat_ctx->temp_solution_D,pcis->vec1_D);CHKERRQ(ierr); 1515 ierr = MatMult(pcis->A_BI,pcis->vec1_D,pcis->vec1_B);CHKERRQ(ierr); 1516 ierr = VecAXPY(mat_ctx->temp_solution_B,-1.0,pcis->vec1_B);CHKERRQ(ierr); 1517 ierr = VecSet(standard_rhs,0.0);CHKERRQ(ierr); 1518 ierr = VecScatterBegin(pcis->global_to_B,mat_ctx->temp_solution_B,standard_rhs,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1519 ierr = VecScatterEnd(pcis->global_to_B,mat_ctx->temp_solution_B,standard_rhs,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1520 /* ierr = PCBDDCScalingRestriction(mat_ctx->pc,standard_rhs,mat_ctx->temp_solution_B);CHKERRQ(ierr); */ 1521 ierr = VecScatterBegin(pcis->global_to_B,standard_rhs,mat_ctx->temp_solution_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1522 ierr = VecScatterEnd(pcis->global_to_B,standard_rhs,mat_ctx->temp_solution_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1523 ierr = VecPointwiseMult(mat_ctx->temp_solution_B,pcis->D,mat_ctx->temp_solution_B);CHKERRQ(ierr); 1524 } 1525 /* BDDC rhs */ 1526 ierr = VecCopy(mat_ctx->temp_solution_B,pcis->vec1_B);CHKERRQ(ierr); 1527 if (pcbddc->switch_static) { 1528 ierr = VecCopy(mat_ctx->temp_solution_D,pcis->vec1_D);CHKERRQ(ierr); 1529 } 1530 /* apply BDDC */ 1531 ierr = PCBDDCApplyInterfacePreconditioner(mat_ctx->pc,PETSC_FALSE);CHKERRQ(ierr); 1532 /* Application of B_delta and assembling of rhs for fetidp fluxes */ 1533 ierr = VecSet(fetidp_flux_rhs,0.0);CHKERRQ(ierr); 1534 ierr = MatMult(mat_ctx->B_delta,pcis->vec1_B,mat_ctx->lambda_local);CHKERRQ(ierr); 1535 ierr = VecScatterBegin(mat_ctx->l2g_lambda,mat_ctx->lambda_local,fetidp_flux_rhs,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1536 ierr = VecScatterEnd (mat_ctx->l2g_lambda,mat_ctx->lambda_local,fetidp_flux_rhs,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1537 /* restore original rhs */ 1538 ierr = VecCopy(pcbddc->original_rhs,standard_rhs);CHKERRQ(ierr); 1539 PetscFunctionReturn(0); 1540 } 1541 1542 #undef __FUNCT__ 1543 #define __FUNCT__ "PCBDDCMatFETIDPGetRHS" 1544 /*@ 1545 PCBDDCMatFETIDPGetRHS - Compute the right-hand side for FETIDP linear system 1546 1547 Collective 1548 1549 Input Parameters: 1550 + fetidp_mat - the FETIDP matrix object obtained by calling PCBDDCCreateFETIDPOperators 1551 . standard_rhs - the right-hand side for your linear system 1552 1553 Output Parameters: 1554 - fetidp_flux_rhs - the right-hand side for the FETIDP linear system 1555 1556 Level: developer 1557 1558 Notes: 1559 1560 .seealso: PCBDDC,PCBDDCCreateFETIDPOperators 1561 @*/ 1562 PetscErrorCode PCBDDCMatFETIDPGetRHS(Mat fetidp_mat, Vec standard_rhs, Vec fetidp_flux_rhs) 1563 { 1564 FETIDPMat_ctx mat_ctx; 1565 PetscErrorCode ierr; 1566 1567 PetscFunctionBegin; 1568 ierr = MatShellGetContext(fetidp_mat,&mat_ctx);CHKERRQ(ierr); 1569 ierr = PetscTryMethod(mat_ctx->pc,"PCBDDCMatFETIDPGetRHS_C",(Mat,Vec,Vec),(fetidp_mat,standard_rhs,fetidp_flux_rhs));CHKERRQ(ierr); 1570 PetscFunctionReturn(0); 1571 } 1572 /* -------------------------------------------------------------------------- */ 1573 1574 #undef __FUNCT__ 1575 #define __FUNCT__ "PCBDDCMatFETIDPGetSolution_BDDC" 1576 static PetscErrorCode PCBDDCMatFETIDPGetSolution_BDDC(Mat fetidp_mat, Vec fetidp_flux_sol, Vec standard_sol) 1577 { 1578 FETIDPMat_ctx mat_ctx; 1579 PC_IS* pcis; 1580 PC_BDDC* pcbddc; 1581 PetscErrorCode ierr; 1582 1583 PetscFunctionBegin; 1584 ierr = MatShellGetContext(fetidp_mat,&mat_ctx);CHKERRQ(ierr); 1585 pcis = (PC_IS*)mat_ctx->pc->data; 1586 pcbddc = (PC_BDDC*)mat_ctx->pc->data; 1587 1588 /* apply B_delta^T */ 1589 ierr = VecScatterBegin(mat_ctx->l2g_lambda,fetidp_flux_sol,mat_ctx->lambda_local,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1590 ierr = VecScatterEnd (mat_ctx->l2g_lambda,fetidp_flux_sol,mat_ctx->lambda_local,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1591 ierr = MatMultTranspose(mat_ctx->B_delta,mat_ctx->lambda_local,pcis->vec1_B);CHKERRQ(ierr); 1592 /* compute rhs for BDDC application */ 1593 ierr = VecAYPX(pcis->vec1_B,-1.0,mat_ctx->temp_solution_B);CHKERRQ(ierr); 1594 if (pcbddc->switch_static) { 1595 ierr = VecCopy(mat_ctx->temp_solution_D,pcis->vec1_D);CHKERRQ(ierr); 1596 } 1597 /* apply BDDC */ 1598 ierr = PCBDDCApplyInterfacePreconditioner(mat_ctx->pc,PETSC_FALSE);CHKERRQ(ierr); 1599 /* put values into standard global vector */ 1600 ierr = VecScatterBegin(pcis->global_to_B,pcis->vec1_B,standard_sol,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1601 ierr = VecScatterEnd (pcis->global_to_B,pcis->vec1_B,standard_sol,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1602 if (!pcbddc->switch_static) { 1603 /* compute values into the interior if solved for the partially subassembled Schur complement */ 1604 ierr = MatMult(pcis->A_IB,pcis->vec1_B,pcis->vec1_D);CHKERRQ(ierr); 1605 ierr = VecAXPY(mat_ctx->temp_solution_D,-1.0,pcis->vec1_D);CHKERRQ(ierr); 1606 ierr = KSPSolve(pcbddc->ksp_D,mat_ctx->temp_solution_D,pcis->vec1_D);CHKERRQ(ierr); 1607 } 1608 ierr = VecScatterBegin(pcis->global_to_D,pcis->vec1_D,standard_sol,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1609 ierr = VecScatterEnd (pcis->global_to_D,pcis->vec1_D,standard_sol,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1610 /* final change of basis if needed 1611 Is also sums the dirichlet part removed during RHS assembling */ 1612 ierr = PCPostSolve_BDDC(mat_ctx->pc,NULL,NULL,standard_sol);CHKERRQ(ierr); 1613 PetscFunctionReturn(0); 1614 } 1615 1616 #undef __FUNCT__ 1617 #define __FUNCT__ "PCBDDCMatFETIDPGetSolution" 1618 /*@ 1619 PCBDDCMatFETIDPGetSolution - Compute the physical solution from the solution of the FETIDP linear system 1620 1621 Collective 1622 1623 Input Parameters: 1624 + fetidp_mat - the FETIDP matrix obtained by calling PCBDDCCreateFETIDPOperators 1625 . fetidp_flux_sol - the solution of the FETIDP linear system 1626 1627 Output Parameters: 1628 - standard_sol - the solution defined on the physical domain 1629 1630 Level: developer 1631 1632 Notes: 1633 1634 .seealso: PCBDDC,PCBDDCCreateFETIDPOperators 1635 @*/ 1636 PetscErrorCode PCBDDCMatFETIDPGetSolution(Mat fetidp_mat, Vec fetidp_flux_sol, Vec standard_sol) 1637 { 1638 FETIDPMat_ctx mat_ctx; 1639 PetscErrorCode ierr; 1640 1641 PetscFunctionBegin; 1642 ierr = MatShellGetContext(fetidp_mat,&mat_ctx);CHKERRQ(ierr); 1643 ierr = PetscTryMethod(mat_ctx->pc,"PCBDDCMatFETIDPGetSolution_C",(Mat,Vec,Vec),(fetidp_mat,fetidp_flux_sol,standard_sol));CHKERRQ(ierr); 1644 PetscFunctionReturn(0); 1645 } 1646 /* -------------------------------------------------------------------------- */ 1647 1648 extern PetscErrorCode FETIDPMatMult(Mat,Vec,Vec); 1649 extern PetscErrorCode FETIDPMatMultTranspose(Mat,Vec,Vec); 1650 extern PetscErrorCode PCBDDCDestroyFETIDPMat(Mat); 1651 extern PetscErrorCode FETIDPPCApply(PC,Vec,Vec); 1652 extern PetscErrorCode FETIDPPCApplyTranspose(PC,Vec,Vec); 1653 extern PetscErrorCode PCBDDCDestroyFETIDPPC(PC); 1654 1655 #undef __FUNCT__ 1656 #define __FUNCT__ "PCBDDCCreateFETIDPOperators_BDDC" 1657 static PetscErrorCode PCBDDCCreateFETIDPOperators_BDDC(PC pc, Mat *fetidp_mat, PC *fetidp_pc) 1658 { 1659 1660 FETIDPMat_ctx fetidpmat_ctx; 1661 Mat newmat; 1662 FETIDPPC_ctx fetidppc_ctx; 1663 PC newpc; 1664 MPI_Comm comm; 1665 PetscErrorCode ierr; 1666 1667 PetscFunctionBegin; 1668 ierr = PetscObjectGetComm((PetscObject)pc,&comm);CHKERRQ(ierr); 1669 /* FETIDP linear matrix */ 1670 ierr = PCBDDCCreateFETIDPMatContext(pc,&fetidpmat_ctx);CHKERRQ(ierr); 1671 ierr = PCBDDCSetupFETIDPMatContext(fetidpmat_ctx);CHKERRQ(ierr); 1672 ierr = MatCreateShell(comm,PETSC_DECIDE,PETSC_DECIDE,fetidpmat_ctx->n_lambda,fetidpmat_ctx->n_lambda,fetidpmat_ctx,&newmat);CHKERRQ(ierr); 1673 ierr = MatShellSetOperation(newmat,MATOP_MULT,(void (*)(void))FETIDPMatMult);CHKERRQ(ierr); 1674 ierr = MatShellSetOperation(newmat,MATOP_MULT_TRANSPOSE,(void (*)(void))FETIDPMatMultTranspose);CHKERRQ(ierr); 1675 ierr = MatShellSetOperation(newmat,MATOP_DESTROY,(void (*)(void))PCBDDCDestroyFETIDPMat);CHKERRQ(ierr); 1676 ierr = MatSetUp(newmat);CHKERRQ(ierr); 1677 /* FETIDP preconditioner */ 1678 ierr = PCBDDCCreateFETIDPPCContext(pc,&fetidppc_ctx);CHKERRQ(ierr); 1679 ierr = PCBDDCSetupFETIDPPCContext(newmat,fetidppc_ctx);CHKERRQ(ierr); 1680 ierr = PCCreate(comm,&newpc);CHKERRQ(ierr); 1681 ierr = PCSetType(newpc,PCSHELL);CHKERRQ(ierr); 1682 ierr = PCShellSetContext(newpc,fetidppc_ctx);CHKERRQ(ierr); 1683 ierr = PCShellSetApply(newpc,FETIDPPCApply);CHKERRQ(ierr); 1684 ierr = PCShellSetApplyTranspose(newpc,FETIDPPCApplyTranspose);CHKERRQ(ierr); 1685 ierr = PCShellSetDestroy(newpc,PCBDDCDestroyFETIDPPC);CHKERRQ(ierr); 1686 ierr = PCSetOperators(newpc,newmat,newmat);CHKERRQ(ierr); 1687 ierr = PCSetUp(newpc);CHKERRQ(ierr); 1688 /* return pointers for objects created */ 1689 *fetidp_mat=newmat; 1690 *fetidp_pc=newpc; 1691 PetscFunctionReturn(0); 1692 } 1693 1694 #undef __FUNCT__ 1695 #define __FUNCT__ "PCBDDCCreateFETIDPOperators" 1696 /*@ 1697 PCBDDCCreateFETIDPOperators - Create operators for FETIDP 1698 1699 Collective 1700 1701 Input Parameters: 1702 + pc - the BDDC preconditioning context already setup 1703 1704 Output Parameters: 1705 . fetidp_mat - shell FETIDP matrix object 1706 . fetidp_pc - shell Dirichlet preconditioner for FETIDP matrix 1707 1708 Options Database Keys: 1709 - -fetidp_fullyredundant: use or not a fully redundant set of Lagrange multipliers 1710 1711 Level: developer 1712 1713 Notes: 1714 Currently the only operation provided for FETIDP matrix is MatMult 1715 1716 .seealso: PCBDDC 1717 @*/ 1718 PetscErrorCode PCBDDCCreateFETIDPOperators(PC pc, Mat *fetidp_mat, PC *fetidp_pc) 1719 { 1720 PetscErrorCode ierr; 1721 1722 PetscFunctionBegin; 1723 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1724 if (pc->setupcalled) { 1725 ierr = PetscUseMethod(pc,"PCBDDCCreateFETIDPOperators_C",(PC,Mat*,PC*),(pc,fetidp_mat,fetidp_pc));CHKERRQ(ierr); 1726 } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"You must call PCSetup_BDDC() first \n"); 1727 PetscFunctionReturn(0); 1728 } 1729 /* -------------------------------------------------------------------------- */ 1730 /*MC 1731 PCBDDC - Balancing Domain Decomposition by Constraints. 1732 1733 An implementation of the BDDC preconditioner based on 1734 1735 .vb 1736 [1] C. R. Dohrmann. "An approximate BDDC preconditioner", Numerical Linear Algebra with Applications Volume 14, Issue 2, pages 149-168, March 2007 1737 [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 1738 [3] J. Mandel, B. Sousedik, C. R. Dohrmann. "Multispace and Multilevel BDDC", http://arxiv.org/abs/0712.3977 1739 .ve 1740 1741 The matrix to be preconditioned (Pmat) must be of type MATIS. 1742 1743 Currently works with MATIS matrices with local Neumann matrices of type MATSEQAIJ, MATSEQBAIJ or MATSEQSBAIJ, either with real or complex numbers. 1744 1745 It also works with unsymmetric and indefinite problems. 1746 1747 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. 1748 1749 Approximate local solvers are automatically adapted for singular linear problems (see [1]) if the user has provided the nullspace using PCBDDCSetNullSpace 1750 1751 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() 1752 1753 Constraints can be customized by attaching a MatNullSpace object to the MATIS matrix via MatSetNearNullSpace(). 1754 1755 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. 1756 1757 The PETSc implementation also supports multilevel BDDC [3]. Coarse grids are partitioned using MatPartitioning object. 1758 1759 Options Database Keys: 1760 1761 . -pc_bddc_use_vertices <1> - use or not vertices in primal space 1762 . -pc_bddc_use_edges <1> - use or not edges in primal space 1763 . -pc_bddc_use_faces <0> - use or not faces in primal space 1764 . -pc_bddc_use_change_of_basis <0> - use change of basis approach (on edges only) 1765 . -pc_bddc_use_change_on_faces <0> - use change of basis approach on faces if change of basis has been requested 1766 . -pc_bddc_switch_static <0> - switches from M_2 to M_3 operator (see reference article [1]) 1767 . -pc_bddc_levels <0> - maximum number of levels for multilevel 1768 . -pc_bddc_coarsening_ratio - H/h ratio at the coarser level 1769 - -pc_bddc_check_level <0> - set verbosity level of debugging output 1770 1771 Options for Dirichlet, Neumann or coarse solver can be set with 1772 .vb 1773 -pc_bddc_dirichlet_ 1774 -pc_bddc_neumann_ 1775 -pc_bddc_coarse_ 1776 .ve 1777 e.g -pc_bddc_dirichlet_ksp_type richardson -pc_bddc_dirichlet_pc_type gamg 1778 1779 When using a multilevel approach, solvers' options at the N-th level can be specified as 1780 .vb 1781 -pc_bddc_dirichlet_lN_ 1782 -pc_bddc_neumann_lN_ 1783 -pc_bddc_coarse_lN_ 1784 .ve 1785 Note that level number ranges from the finest 0 to the coarsest N. 1786 1787 Level: intermediate 1788 1789 Developer notes: 1790 1791 New deluxe scaling operator will be available soon. 1792 1793 Contributed by Stefano Zampini 1794 1795 .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PC, MATIS 1796 M*/ 1797 1798 #undef __FUNCT__ 1799 #define __FUNCT__ "PCCreate_BDDC" 1800 PETSC_EXTERN PetscErrorCode PCCreate_BDDC(PC pc) 1801 { 1802 PetscErrorCode ierr; 1803 PC_BDDC *pcbddc; 1804 1805 PetscFunctionBegin; 1806 /* Creates the private data structure for this preconditioner and attach it to the PC object. */ 1807 ierr = PetscNewLog(pc,&pcbddc);CHKERRQ(ierr); 1808 pc->data = (void*)pcbddc; 1809 1810 /* create PCIS data structure */ 1811 ierr = PCISCreate(pc);CHKERRQ(ierr); 1812 1813 /* BDDC customization */ 1814 pcbddc->use_local_adj = PETSC_TRUE; 1815 pcbddc->use_vertices = PETSC_TRUE; 1816 pcbddc->use_edges = PETSC_TRUE; 1817 pcbddc->use_faces = PETSC_FALSE; 1818 pcbddc->use_change_of_basis = PETSC_FALSE; 1819 pcbddc->use_change_on_faces = PETSC_FALSE; 1820 pcbddc->switch_static = PETSC_FALSE; 1821 pcbddc->use_nnsp_true = PETSC_FALSE; /* not yet exposed */ 1822 pcbddc->dbg_flag = 0; 1823 /* private */ 1824 pcbddc->issym = PETSC_FALSE; 1825 pcbddc->BtoNmap = 0; 1826 pcbddc->local_primal_size = 0; 1827 pcbddc->n_vertices = 0; 1828 pcbddc->n_actual_vertices = 0; 1829 pcbddc->n_constraints = 0; 1830 pcbddc->primal_indices_local_idxs = 0; 1831 pcbddc->recompute_topography = PETSC_FALSE; 1832 pcbddc->coarse_size = -1; 1833 pcbddc->new_primal_space = PETSC_FALSE; 1834 pcbddc->new_primal_space_local = PETSC_FALSE; 1835 pcbddc->global_primal_indices = 0; 1836 pcbddc->onearnullspace = 0; 1837 pcbddc->onearnullvecs_state = 0; 1838 pcbddc->user_primal_vertices = 0; 1839 pcbddc->NullSpace = 0; 1840 pcbddc->temp_solution = 0; 1841 pcbddc->original_rhs = 0; 1842 pcbddc->local_mat = 0; 1843 pcbddc->ChangeOfBasisMatrix = 0; 1844 pcbddc->user_ChangeOfBasisMatrix = 0; 1845 pcbddc->new_global_mat = 0; 1846 pcbddc->coarse_vec = 0; 1847 pcbddc->coarse_rhs = 0; 1848 pcbddc->coarse_ksp = 0; 1849 pcbddc->coarse_phi_B = 0; 1850 pcbddc->coarse_phi_D = 0; 1851 pcbddc->coarse_psi_B = 0; 1852 pcbddc->coarse_psi_D = 0; 1853 pcbddc->vec1_P = 0; 1854 pcbddc->vec1_R = 0; 1855 pcbddc->vec2_R = 0; 1856 pcbddc->local_auxmat1 = 0; 1857 pcbddc->local_auxmat2 = 0; 1858 pcbddc->R_to_B = 0; 1859 pcbddc->R_to_D = 0; 1860 pcbddc->ksp_D = 0; 1861 pcbddc->ksp_R = 0; 1862 pcbddc->NeumannBoundaries = 0; 1863 pcbddc->NeumannBoundariesLocal = 0; 1864 pcbddc->DirichletBoundaries = 0; 1865 pcbddc->DirichletBoundariesLocal = 0; 1866 pcbddc->user_provided_isfordofs = PETSC_FALSE; 1867 pcbddc->n_ISForDofs = 0; 1868 pcbddc->n_ISForDofsLocal = 0; 1869 pcbddc->ISForDofs = 0; 1870 pcbddc->ISForDofsLocal = 0; 1871 pcbddc->ConstraintMatrix = 0; 1872 pcbddc->use_exact_dirichlet_trick = PETSC_TRUE; 1873 pcbddc->coarse_loc_to_glob = 0; 1874 pcbddc->coarsening_ratio = 8; 1875 pcbddc->current_level = 0; 1876 pcbddc->max_levels = 0; 1877 pcbddc->use_coarse_estimates = PETSC_FALSE; 1878 pcbddc->coarse_subassembling = 0; 1879 pcbddc->coarse_subassembling_init = 0; 1880 1881 /* create local graph structure */ 1882 ierr = PCBDDCGraphCreate(&pcbddc->mat_graph);CHKERRQ(ierr); 1883 1884 /* scaling */ 1885 pcbddc->work_scaling = 0; 1886 pcbddc->use_deluxe_scaling = PETSC_FALSE; 1887 pcbddc->deluxe_threshold = 1; 1888 pcbddc->deluxe_rebuild = PETSC_FALSE; 1889 pcbddc->deluxe_layers = -1; 1890 pcbddc->deluxe_use_useradj = PETSC_FALSE; 1891 pcbddc->deluxe_compute_rowadj = PETSC_TRUE; 1892 1893 /* function pointers */ 1894 pc->ops->apply = PCApply_BDDC; 1895 pc->ops->applytranspose = PCApplyTranspose_BDDC; 1896 pc->ops->setup = PCSetUp_BDDC; 1897 pc->ops->destroy = PCDestroy_BDDC; 1898 pc->ops->setfromoptions = PCSetFromOptions_BDDC; 1899 pc->ops->view = 0; 1900 pc->ops->applyrichardson = 0; 1901 pc->ops->applysymmetricleft = 0; 1902 pc->ops->applysymmetricright = 0; 1903 pc->ops->presolve = PCPreSolve_BDDC; 1904 pc->ops->postsolve = PCPostSolve_BDDC; 1905 1906 /* composing function */ 1907 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetChangeOfBasisMat_C",PCBDDCSetChangeOfBasisMat_BDDC);CHKERRQ(ierr); 1908 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetPrimalVerticesLocalIS_C",PCBDDCSetPrimalVerticesLocalIS_BDDC);CHKERRQ(ierr); 1909 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetCoarseningRatio_C",PCBDDCSetCoarseningRatio_BDDC);CHKERRQ(ierr); 1910 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetLevel_C",PCBDDCSetLevel_BDDC);CHKERRQ(ierr); 1911 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetUseExactDirichlet_C",PCBDDCSetUseExactDirichlet_BDDC);CHKERRQ(ierr); 1912 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetLevels_C",PCBDDCSetLevels_BDDC);CHKERRQ(ierr); 1913 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetNullSpace_C",PCBDDCSetNullSpace_BDDC);CHKERRQ(ierr); 1914 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetDirichletBoundaries_C",PCBDDCSetDirichletBoundaries_BDDC);CHKERRQ(ierr); 1915 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetDirichletBoundariesLocal_C",PCBDDCSetDirichletBoundariesLocal_BDDC);CHKERRQ(ierr); 1916 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetNeumannBoundaries_C",PCBDDCSetNeumannBoundaries_BDDC);CHKERRQ(ierr); 1917 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetNeumannBoundariesLocal_C",PCBDDCSetNeumannBoundariesLocal_BDDC);CHKERRQ(ierr); 1918 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCGetDirichletBoundaries_C",PCBDDCGetDirichletBoundaries_BDDC);CHKERRQ(ierr); 1919 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCGetDirichletBoundariesLocal_C",PCBDDCGetDirichletBoundariesLocal_BDDC);CHKERRQ(ierr); 1920 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCGetNeumannBoundaries_C",PCBDDCGetNeumannBoundaries_BDDC);CHKERRQ(ierr); 1921 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCGetNeumannBoundariesLocal_C",PCBDDCGetNeumannBoundariesLocal_BDDC);CHKERRQ(ierr); 1922 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetDofsSplitting_C",PCBDDCSetDofsSplitting_BDDC);CHKERRQ(ierr); 1923 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetDofsSplittingLocal_C",PCBDDCSetDofsSplittingLocal_BDDC);CHKERRQ(ierr); 1924 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetLocalAdjacencyGraph_C",PCBDDCSetLocalAdjacencyGraph_BDDC);CHKERRQ(ierr); 1925 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCCreateFETIDPOperators_C",PCBDDCCreateFETIDPOperators_BDDC);CHKERRQ(ierr); 1926 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCMatFETIDPGetRHS_C",PCBDDCMatFETIDPGetRHS_BDDC);CHKERRQ(ierr); 1927 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCMatFETIDPGetSolution_C",PCBDDCMatFETIDPGetSolution_BDDC);CHKERRQ(ierr); 1928 PetscFunctionReturn(0); 1929 } 1930 1931