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