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