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