1 #include <petsc/private/kspimpl.h> /*I <petscksp.h> I*/ 2 #include <../src/ksp/pc/impls/bddc/bddc.h> 3 #include <../src/ksp/pc/impls/bddc/bddcprivate.h> 4 5 /* 6 This file implements the FETI-DP method in PETSc as part of KSP. 7 */ 8 typedef struct { 9 KSP parentksp; 10 } KSP_FETIDPMon; 11 12 typedef struct { 13 KSP innerksp; /* the KSP for the Lagrange multipliers */ 14 PC innerbddc; /* the inner BDDC object */ 15 PetscBool fully_redundant; /* true for using a fully redundant set of multipliers */ 16 PetscBool userbddc; /* true if the user provided the PCBDDC object */ 17 PetscBool saddlepoint; /* support for saddle point problems */ 18 IS pP; /* index set for interface pressure variables */ 19 Vec rhs_flip; /* see KSPFETIDPSetUpOperators */ 20 KSP_FETIDPMon *monctx; /* monitor context, used to pass user defined monitors 21 in the physical space */ 22 PetscObjectState matstate; /* these are needed just in the saddle point case */ 23 PetscObjectState matnnzstate; /* where we are going to use MatZeroRows on pmat */ 24 PetscBool statechanged; 25 PetscBool check; 26 } KSP_FETIDP; 27 28 #undef __FUNCT__ 29 #define __FUNCT__ "KSPFETIDPSetPressureOperators_FETIDP" 30 static PetscErrorCode KSPFETIDPSetPressureOperators_FETIDP(KSP ksp, Mat A, Mat P) 31 { 32 KSP_FETIDP *fetidp = (KSP_FETIDP*)ksp->data; 33 PetscErrorCode ierr; 34 35 PetscFunctionBegin; 36 if (A || P) fetidp->saddlepoint = PETSC_TRUE; 37 ierr = PetscObjectCompose((PetscObject)fetidp->innerbddc,"__KSPFETIDP_PAmat",(PetscObject)A);CHKERRQ(ierr); 38 ierr = PetscObjectCompose((PetscObject)fetidp->innerbddc,"__KSPFETIDP_PPmat",(PetscObject)P);CHKERRQ(ierr); 39 PetscFunctionReturn(0); 40 } 41 42 #undef __FUNCT__ 43 #define __FUNCT__ "KSPFETIDPSetPressureOperators" 44 /*@ 45 KSPFETIDPSetPressureOperators - Sets the operators used to setup the pressure preconditioner for saddle point FETI-DP. 46 47 Collective on KSP 48 49 Input Parameters: 50 + ksp - the FETI-DP Krylov solver 51 . A - the linear operator defining the pressure problem 52 - P - the linear operator to be preconditioned 53 54 Level: advanced 55 56 Notes: The operators can be either passed in monolithic global ordering or in interface pressure ordering. 57 In the latter case, the interface pressure ordering of dofs needs to satisfy 58 pid_1 < pid_2 iff gid_1 < gid_2 59 where pid_1 and pid_2 are two different pressure dof numbers and gid_1 and gid_2 the corresponding 60 id in the global ordering. 61 62 .seealso: MATIS, PCBDDC, KSPFETIDPGetInnerBDDC, KSPFETIDPGetInnerKSP, KSPSetOperators 63 @*/ 64 PetscErrorCode KSPFETIDPSetPressureOperators(KSP ksp, Mat A, Mat P) 65 { 66 PetscErrorCode ierr; 67 68 PetscFunctionBegin; 69 PetscValidHeaderSpecific(ksp,KSP_CLASSID,1); 70 if (A) PetscValidHeaderSpecific(A,MAT_CLASSID,2); 71 if (P) PetscValidHeaderSpecific(P,MAT_CLASSID,3); 72 ierr = PetscTryMethod(ksp,"KSPFETIDPSetPressureOperators_C",(KSP,Mat,Mat),(ksp,A,P));CHKERRQ(ierr); 73 PetscFunctionReturn(0); 74 } 75 76 #undef __FUNCT__ 77 #define __FUNCT__ "KSPFETIDPGetInnerKSP_FETIDP" 78 static PetscErrorCode KSPFETIDPGetInnerKSP_FETIDP(KSP ksp, KSP* innerksp) 79 { 80 KSP_FETIDP *fetidp = (KSP_FETIDP*)ksp->data; 81 82 PetscFunctionBegin; 83 *innerksp = fetidp->innerksp; 84 PetscFunctionReturn(0); 85 } 86 87 #undef __FUNCT__ 88 #define __FUNCT__ "KSPFETIDPGetInnerKSP" 89 /*@ 90 KSPFETIDPGetInnerKSP - Gets the KSP object for the Lagrange multipliers 91 92 Input Parameters: 93 + ksp - the FETI-DP KSP 94 - innerksp - the KSP for the multipliers 95 96 Level: advanced 97 98 Notes: 99 100 .seealso: MATIS, PCBDDC, KSPFETIDPSetInnerBDDC, KSPFETIDPGetInnerBDDC 101 @*/ 102 PetscErrorCode KSPFETIDPGetInnerKSP(KSP ksp, KSP* innerksp) 103 { 104 PetscErrorCode ierr; 105 106 PetscFunctionBegin; 107 PetscValidHeaderSpecific(ksp,KSP_CLASSID,1); 108 PetscValidPointer(innerksp,2); 109 ierr = PetscUseMethod(ksp,"KSPFETIDPGetInnerKSP_C",(KSP,KSP*),(ksp,innerksp));CHKERRQ(ierr); 110 PetscFunctionReturn(0); 111 } 112 113 #undef __FUNCT__ 114 #define __FUNCT__ "KSPFETIDPGetInnerBDDC_FETIDP" 115 static PetscErrorCode KSPFETIDPGetInnerBDDC_FETIDP(KSP ksp, PC* pc) 116 { 117 KSP_FETIDP *fetidp = (KSP_FETIDP*)ksp->data; 118 119 PetscFunctionBegin; 120 *pc = fetidp->innerbddc; 121 PetscFunctionReturn(0); 122 } 123 124 #undef __FUNCT__ 125 #define __FUNCT__ "KSPFETIDPGetInnerBDDC" 126 /*@ 127 KSPFETIDPGetInnerBDDC - Gets the BDDC preconditioner used to setup the FETI-DP matrix for the Lagrange multipliers 128 129 Input Parameters: 130 + ksp - the FETI-DP Krylov solver 131 - pc - the BDDC preconditioner 132 133 Level: advanced 134 135 Notes: 136 137 .seealso: MATIS, PCBDDC, KSPFETIDPSetInnerBDDC, KSPFETIDPGetInnerKSP 138 @*/ 139 PetscErrorCode KSPFETIDPGetInnerBDDC(KSP ksp, PC* pc) 140 { 141 PetscErrorCode ierr; 142 143 PetscFunctionBegin; 144 PetscValidHeaderSpecific(ksp,KSP_CLASSID,1); 145 PetscValidPointer(pc,2); 146 ierr = PetscUseMethod(ksp,"KSPFETIDPGetInnerBDDC_C",(KSP,PC*),(ksp,pc));CHKERRQ(ierr); 147 PetscFunctionReturn(0); 148 } 149 150 #undef __FUNCT__ 151 #define __FUNCT__ "KSPFETIDPSetInnerBDDC_FETIDP" 152 static PetscErrorCode KSPFETIDPSetInnerBDDC_FETIDP(KSP ksp, PC pc) 153 { 154 KSP_FETIDP *fetidp = (KSP_FETIDP*)ksp->data; 155 PetscErrorCode ierr; 156 157 PetscFunctionBegin; 158 ierr = PetscObjectReference((PetscObject)pc);CHKERRQ(ierr); 159 ierr = PCDestroy(&fetidp->innerbddc);CHKERRQ(ierr); 160 fetidp->innerbddc = pc; 161 fetidp->userbddc = PETSC_TRUE; 162 PetscFunctionReturn(0); 163 } 164 165 #undef __FUNCT__ 166 #define __FUNCT__ "KSPFETIDPSetInnerBDDC" 167 /*@ 168 KSPFETIDPSetInnerBDDC - Sets the BDDC preconditioner used to setup the FETI-DP matrix for the Lagrange multipliers 169 170 Collective on KSP 171 172 Input Parameters: 173 + ksp - the FETI-DP Krylov solver 174 - pc - the BDDC preconditioner 175 176 Level: advanced 177 178 Notes: 179 180 .seealso: MATIS, PCBDDC, KSPFETIDPGetInnerBDDC, KSPFETIDPGetInnerKSP 181 @*/ 182 PetscErrorCode KSPFETIDPSetInnerBDDC(KSP ksp, PC pc) 183 { 184 PetscBool isbddc; 185 PetscErrorCode ierr; 186 187 PetscFunctionBegin; 188 PetscValidHeaderSpecific(ksp,KSP_CLASSID,1); 189 PetscValidHeaderSpecific(pc,PC_CLASSID,2); 190 ierr = PetscObjectTypeCompare((PetscObject)pc,PCBDDC,&isbddc);CHKERRQ(ierr); 191 if (!isbddc) SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_ARG_WRONG,"KSPFETIDPSetInnerBDDC need a PCBDDC preconditioner"); 192 ierr = PetscTryMethod(ksp,"KSPFETIDPSetInnerBDDC_C",(KSP,PC),(ksp,pc));CHKERRQ(ierr); 193 PetscFunctionReturn(0); 194 } 195 196 #undef __FUNCT__ 197 #define __FUNCT__ "KSPBuildSolution_FETIDP" 198 static PetscErrorCode KSPBuildSolution_FETIDP(KSP ksp,Vec v,Vec *V) 199 { 200 KSP_FETIDP *fetidp = (KSP_FETIDP*)ksp->data; 201 Mat F; 202 Vec Xl; 203 PetscErrorCode ierr; 204 205 PetscFunctionBegin; 206 ierr = KSPGetOperators(fetidp->innerksp,&F,NULL);CHKERRQ(ierr); 207 ierr = KSPBuildSolution(fetidp->innerksp,NULL,&Xl);CHKERRQ(ierr); 208 if (v) { 209 ierr = PCBDDCMatFETIDPGetSolution(F,Xl,v);CHKERRQ(ierr); 210 *V = v; 211 } else { 212 ierr = PCBDDCMatFETIDPGetSolution(F,Xl,*V);CHKERRQ(ierr); 213 } 214 PetscFunctionReturn(0); 215 } 216 217 #undef __FUNCT__ 218 #define __FUNCT__ "KSPMonitor_FETIDP" 219 static PetscErrorCode KSPMonitor_FETIDP(KSP ksp,PetscInt it,PetscReal rnorm,void* ctx) 220 { 221 KSP_FETIDPMon *monctx = (KSP_FETIDPMon*)ctx; 222 PetscErrorCode ierr; 223 224 PetscFunctionBegin; 225 ierr = KSPMonitor(monctx->parentksp,it,rnorm);CHKERRQ(ierr); 226 PetscFunctionReturn(0); 227 } 228 229 #undef __FUNCT__ 230 #define __FUNCT__ "KSPComputeEigenvalues_FETIDP" 231 static PetscErrorCode KSPComputeEigenvalues_FETIDP(KSP ksp,PetscInt nmax,PetscReal *r,PetscReal *c,PetscInt *neig) 232 { 233 KSP_FETIDP *fetidp = (KSP_FETIDP*)ksp->data; 234 PetscErrorCode ierr; 235 236 PetscFunctionBegin; 237 ierr = KSPComputeEigenvalues(fetidp->innerksp,nmax,r,c,neig);CHKERRQ(ierr); 238 PetscFunctionReturn(0); 239 } 240 241 #undef __FUNCT__ 242 #define __FUNCT__ "KSPComputeExtremeSingularValues_FETIDP" 243 static PetscErrorCode KSPComputeExtremeSingularValues_FETIDP(KSP ksp,PetscReal *emax,PetscReal *emin) 244 { 245 KSP_FETIDP *fetidp = (KSP_FETIDP*)ksp->data; 246 PetscErrorCode ierr; 247 248 PetscFunctionBegin; 249 ierr = KSPComputeExtremeSingularValues(fetidp->innerksp,emax,emin);CHKERRQ(ierr); 250 PetscFunctionReturn(0); 251 } 252 253 #undef __FUNCT__ 254 #define __FUNCT__ "KSPFETIDPCheckOperators" 255 static PetscErrorCode KSPFETIDPCheckOperators(KSP ksp, PetscViewer viewer) 256 { 257 KSP_FETIDP *fetidp = (KSP_FETIDP*)ksp->data; 258 PC_BDDC *pcbddc = (PC_BDDC*)fetidp->innerbddc->data; 259 PC_IS *pcis = (PC_IS*)fetidp->innerbddc->data; 260 Mat_IS *matis = (Mat_IS*)fetidp->innerbddc->pmat->data; 261 Mat F; 262 FETIDPMat_ctx fetidpmat_ctx; 263 Vec test_vec,test_vec_p = NULL,fetidp_global; 264 IS dirdofs,isvert; 265 MPI_Comm comm = PetscObjectComm((PetscObject)ksp); 266 PetscScalar sval,*array; 267 PetscReal val,rval; 268 const PetscInt *vertex_indices; 269 PetscInt i,n_vertices; 270 PetscBool isascii; 271 PetscErrorCode ierr; 272 273 PetscFunctionBegin; 274 PetscCheckSameComm(ksp,1,viewer,2); 275 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);CHKERRQ(ierr); 276 if (!isascii) SETERRQ(comm,PETSC_ERR_SUP,"Unsupported viewer"); 277 ierr = PetscViewerASCIIPrintf(viewer,"----------FETI-DP MAT --------------\n");CHKERRQ(ierr); 278 ierr = PetscViewerASCIIAddTab(viewer,2);CHKERRQ(ierr); 279 ierr = KSPGetOperators(fetidp->innerksp,&F,NULL);CHKERRQ(ierr); 280 ierr = PetscViewerPushFormat(viewer,PETSC_VIEWER_ASCII_INFO);CHKERRQ(ierr); 281 ierr = MatView(F,viewer);CHKERRQ(ierr); 282 ierr = PetscViewerPopFormat(viewer);CHKERRQ(ierr); 283 ierr = PetscViewerASCIISubtractTab(viewer,2);CHKERRQ(ierr); 284 ierr = MatShellGetContext(F,(void**)&fetidpmat_ctx);CHKERRQ(ierr); 285 ierr = PetscViewerASCIIPrintf(viewer,"----------FETI-DP TESTS--------------\n");CHKERRQ(ierr); 286 ierr = PetscViewerASCIIPrintf(viewer,"All tests should return zero!\n");CHKERRQ(ierr); 287 ierr = PetscViewerASCIIPrintf(viewer,"FETIDP MAT context in the ");CHKERRQ(ierr); 288 if (fetidp->fully_redundant) { 289 ierr = PetscViewerASCIIPrintf(viewer,"fully redundant case for lagrange multipliers.\n");CHKERRQ(ierr); 290 } else { 291 ierr = PetscViewerASCIIPrintf(viewer,"Non-fully redundant case for lagrange multiplier.\n");CHKERRQ(ierr); 292 } 293 ierr = PetscViewerFlush(viewer);CHKERRQ(ierr); 294 295 /* Get Vertices used to define the BDDC */ 296 ierr = PCBDDCGraphGetCandidatesIS(pcbddc->mat_graph,NULL,NULL,NULL,NULL,&isvert);CHKERRQ(ierr); 297 ierr = ISGetLocalSize(isvert,&n_vertices);CHKERRQ(ierr); 298 ierr = ISGetIndices(isvert,&vertex_indices);CHKERRQ(ierr); 299 300 /******************************************************************/ 301 /* TEST A/B: Test numbering of global fetidp dofs */ 302 /******************************************************************/ 303 ierr = MatCreateVecs(F,&fetidp_global,NULL);CHKERRQ(ierr); 304 ierr = VecDuplicate(fetidpmat_ctx->lambda_local,&test_vec);CHKERRQ(ierr); 305 ierr = VecSet(fetidp_global,1.0);CHKERRQ(ierr); 306 ierr = VecSet(test_vec,1.);CHKERRQ(ierr); 307 ierr = VecScatterBegin(fetidpmat_ctx->l2g_lambda,fetidp_global,fetidpmat_ctx->lambda_local,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 308 ierr = VecScatterEnd(fetidpmat_ctx->l2g_lambda,fetidp_global,fetidpmat_ctx->lambda_local,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 309 if (fetidpmat_ctx->l2g_p) { 310 ierr = VecDuplicate(fetidpmat_ctx->vP,&test_vec_p);CHKERRQ(ierr); 311 ierr = VecSet(test_vec_p,1.);CHKERRQ(ierr); 312 ierr = VecScatterBegin(fetidpmat_ctx->l2g_p,fetidp_global,fetidpmat_ctx->vP,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 313 ierr = VecScatterEnd(fetidpmat_ctx->l2g_p,fetidp_global,fetidpmat_ctx->vP,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 314 } 315 ierr = VecAXPY(test_vec,-1.0,fetidpmat_ctx->lambda_local);CHKERRQ(ierr); 316 ierr = VecNorm(test_vec,NORM_INFINITY,&val);CHKERRQ(ierr); 317 ierr = VecDestroy(&test_vec);CHKERRQ(ierr); 318 ierr = MPI_Reduce(&val,&rval,1,MPIU_REAL,MPI_MAX,0,comm);CHKERRQ(ierr); 319 ierr = PetscViewerASCIIPrintf(viewer,"A: CHECK glob to loc: % 1.14e\n",rval);CHKERRQ(ierr); 320 321 if (fetidpmat_ctx->l2g_p) { 322 ierr = VecAXPY(test_vec_p,-1.0,fetidpmat_ctx->vP);CHKERRQ(ierr); 323 ierr = VecNorm(test_vec_p,NORM_INFINITY,&val);CHKERRQ(ierr); 324 ierr = MPI_Reduce(&val,&rval,1,MPIU_REAL,MPI_MAX,0,comm);CHKERRQ(ierr); 325 ierr = PetscViewerASCIIPrintf(viewer,"A: CHECK glob to loc (p): % 1.14e\n",rval);CHKERRQ(ierr); 326 } 327 328 if (fetidp->fully_redundant) { 329 ierr = VecSet(fetidp_global,0.0);CHKERRQ(ierr); 330 ierr = VecSet(fetidpmat_ctx->lambda_local,0.5);CHKERRQ(ierr); 331 ierr = VecScatterBegin(fetidpmat_ctx->l2g_lambda,fetidpmat_ctx->lambda_local,fetidp_global,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 332 ierr = VecScatterEnd(fetidpmat_ctx->l2g_lambda,fetidpmat_ctx->lambda_local,fetidp_global,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 333 ierr = VecSum(fetidp_global,&sval);CHKERRQ(ierr); 334 val = PetscRealPart(sval)-fetidpmat_ctx->n_lambda; 335 ierr = MPI_Reduce(&val,&rval,1,MPIU_REAL,MPI_MAX,0,comm);CHKERRQ(ierr); 336 ierr = PetscViewerASCIIPrintf(viewer,"B: CHECK loc to glob: % 1.14e\n",rval);CHKERRQ(ierr); 337 } 338 339 if (fetidpmat_ctx->l2g_p) { 340 ierr = VecSet(pcis->vec1_N,1.0);CHKERRQ(ierr); 341 ierr = VecSet(pcis->vec1_global,0.0);CHKERRQ(ierr); 342 ierr = VecScatterBegin(matis->rctx,pcis->vec1_N,pcis->vec1_global,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 343 ierr = VecScatterEnd(matis->rctx,pcis->vec1_N,pcis->vec1_global,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 344 345 ierr = VecSet(fetidp_global,0.0);CHKERRQ(ierr); 346 ierr = VecSet(fetidpmat_ctx->vP,-1.0);CHKERRQ(ierr); 347 ierr = VecScatterBegin(fetidpmat_ctx->l2g_p,fetidpmat_ctx->vP,fetidp_global,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 348 ierr = VecScatterEnd(fetidpmat_ctx->l2g_p,fetidpmat_ctx->vP,fetidp_global,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 349 ierr = VecScatterBegin(fetidpmat_ctx->g2g_p,fetidp_global,pcis->vec1_global,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 350 ierr = VecScatterEnd(fetidpmat_ctx->g2g_p,fetidp_global,pcis->vec1_global,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 351 ierr = VecScatterBegin(fetidpmat_ctx->g2g_p,pcis->vec1_global,fetidp_global,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 352 ierr = VecScatterEnd(fetidpmat_ctx->g2g_p,pcis->vec1_global,fetidp_global,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 353 ierr = VecSum(fetidp_global,&sval);CHKERRQ(ierr); 354 val = PetscRealPart(sval); 355 ierr = MPI_Reduce(&val,&rval,1,MPIU_REAL,MPI_MAX,0,comm);CHKERRQ(ierr); 356 ierr = PetscViewerASCIIPrintf(viewer,"B: CHECK loc to glob (p): % 1.14e\n",rval);CHKERRQ(ierr); 357 } 358 359 /******************************************************************/ 360 /* TEST C: It should hold B_delta*w=0, w\in\widehat{W} */ 361 /* This is the meaning of the B matrix */ 362 /******************************************************************/ 363 364 ierr = VecSetRandom(pcis->vec1_N,NULL);CHKERRQ(ierr); 365 ierr = VecSet(pcis->vec1_global,0.0);CHKERRQ(ierr); 366 ierr = VecScatterBegin(matis->rctx,pcis->vec1_N,pcis->vec1_global,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 367 ierr = VecScatterEnd(matis->rctx,pcis->vec1_N,pcis->vec1_global,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 368 ierr = VecScatterBegin(matis->rctx,pcis->vec1_global,pcis->vec1_N,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 369 ierr = VecScatterEnd(matis->rctx,pcis->vec1_global,pcis->vec1_N,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 370 ierr = VecScatterBegin(pcis->N_to_B,pcis->vec1_N,pcis->vec1_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 371 ierr = VecScatterEnd(pcis->N_to_B,pcis->vec1_N,pcis->vec1_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 372 /* Action of B_delta */ 373 ierr = MatMult(fetidpmat_ctx->B_delta,pcis->vec1_B,fetidpmat_ctx->lambda_local);CHKERRQ(ierr); 374 ierr = VecSet(fetidp_global,0.0);CHKERRQ(ierr); 375 ierr = VecScatterBegin(fetidpmat_ctx->l2g_lambda,fetidpmat_ctx->lambda_local,fetidp_global,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 376 ierr = VecScatterEnd(fetidpmat_ctx->l2g_lambda,fetidpmat_ctx->lambda_local,fetidp_global,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 377 ierr = VecNorm(fetidp_global,NORM_INFINITY,&val);CHKERRQ(ierr); 378 ierr = PetscViewerASCIIPrintf(viewer,"C: CHECK infty norm of B_delta*w (w continuous): % 1.14e\n",val);CHKERRQ(ierr); 379 380 /******************************************************************/ 381 /* TEST D: It should hold E_Dw = w - P_Dw w\in\widetilde{W} */ 382 /* E_D = R_D^TR */ 383 /* P_D = B_{D,delta}^T B_{delta} */ 384 /* eq.44 Mandel Tezaur and Dohrmann 2005 */ 385 /******************************************************************/ 386 387 /* compute a random vector in \widetilde{W} */ 388 ierr = VecSetRandom(pcis->vec1_N,NULL);CHKERRQ(ierr); 389 /* set zero at vertices and essential dofs */ 390 ierr = VecGetArray(pcis->vec1_N,&array);CHKERRQ(ierr); 391 for (i=0;i<n_vertices;i++) array[vertex_indices[i]] = 0.0; 392 ierr = PCBDDCGraphGetDirichletDofs(pcbddc->mat_graph,&dirdofs);CHKERRQ(ierr); 393 if (dirdofs) { 394 const PetscInt *idxs; 395 PetscInt ndir; 396 397 ierr = ISGetLocalSize(dirdofs,&ndir);CHKERRQ(ierr); 398 ierr = ISGetIndices(dirdofs,&idxs);CHKERRQ(ierr); 399 for (i=0;i<ndir;i++) array[idxs[i]] = 0.0; 400 ierr = ISRestoreIndices(dirdofs,&idxs);CHKERRQ(ierr); 401 } 402 ierr = VecRestoreArray(pcis->vec1_N,&array);CHKERRQ(ierr); 403 /* store w for final comparison */ 404 ierr = VecDuplicate(pcis->vec1_B,&test_vec);CHKERRQ(ierr); 405 ierr = VecScatterBegin(pcis->N_to_B,pcis->vec1_N,test_vec,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 406 ierr = VecScatterEnd(pcis->N_to_B,pcis->vec1_N,test_vec,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 407 408 /* Jump operator P_D : results stored in pcis->vec1_B */ 409 /* Action of B_delta */ 410 ierr = MatMult(fetidpmat_ctx->B_delta,test_vec,fetidpmat_ctx->lambda_local);CHKERRQ(ierr); 411 ierr = VecSet(fetidp_global,0.0);CHKERRQ(ierr); 412 ierr = VecScatterBegin(fetidpmat_ctx->l2g_lambda,fetidpmat_ctx->lambda_local,fetidp_global,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 413 ierr = VecScatterEnd(fetidpmat_ctx->l2g_lambda,fetidpmat_ctx->lambda_local,fetidp_global,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 414 /* Action of B_Ddelta^T */ 415 ierr = VecScatterBegin(fetidpmat_ctx->l2g_lambda,fetidp_global,fetidpmat_ctx->lambda_local,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 416 ierr = VecScatterEnd(fetidpmat_ctx->l2g_lambda,fetidp_global,fetidpmat_ctx->lambda_local,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 417 ierr = MatMultTranspose(fetidpmat_ctx->B_Ddelta,fetidpmat_ctx->lambda_local,pcis->vec1_B);CHKERRQ(ierr); 418 419 /* Average operator E_D : results stored in pcis->vec2_B */ 420 ierr = PCBDDCScalingExtension(fetidpmat_ctx->pc,test_vec,pcis->vec1_global);CHKERRQ(ierr); 421 ierr = VecScatterBegin(pcis->global_to_B,pcis->vec1_global,pcis->vec2_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 422 ierr = VecScatterEnd(pcis->global_to_B,pcis->vec1_global,pcis->vec2_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 423 424 /* test E_D=I-P_D */ 425 ierr = VecAXPY(pcis->vec1_B,1.0,pcis->vec2_B);CHKERRQ(ierr); 426 ierr = VecAXPY(pcis->vec1_B,-1.0,test_vec);CHKERRQ(ierr); 427 ierr = VecNorm(pcis->vec1_B,NORM_INFINITY,&val);CHKERRQ(ierr); 428 ierr = VecDestroy(&test_vec);CHKERRQ(ierr); 429 ierr = MPI_Reduce(&val,&rval,1,MPIU_REAL,MPI_MAX,0,comm);CHKERRQ(ierr); 430 ierr = PetscViewerASCIIPrintf(viewer,"D: CHECK infty norm of E_D + P_D - I: % 1.14e\n",PetscGlobalRank,val);CHKERRQ(ierr); 431 432 /******************************************************************/ 433 /* TEST E: It should hold R_D^TP_Dw=0 w\in\widetilde{W} */ 434 /* eq.48 Mandel Tezaur and Dohrmann 2005 */ 435 /******************************************************************/ 436 437 ierr = VecSetRandom(pcis->vec1_N,NULL);CHKERRQ(ierr); 438 /* set zero at vertices and essential dofs */ 439 ierr = VecGetArray(pcis->vec1_N,&array);CHKERRQ(ierr); 440 for (i=0;i<n_vertices;i++) array[vertex_indices[i]] = 0.0; 441 if (dirdofs) { 442 const PetscInt *idxs; 443 PetscInt ndir; 444 445 ierr = ISGetLocalSize(dirdofs,&ndir);CHKERRQ(ierr); 446 ierr = ISGetIndices(dirdofs,&idxs);CHKERRQ(ierr); 447 for (i=0;i<ndir;i++) array[idxs[i]] = 0.0; 448 ierr = ISRestoreIndices(dirdofs,&idxs);CHKERRQ(ierr); 449 } 450 ierr = VecRestoreArray(pcis->vec1_N,&array);CHKERRQ(ierr); 451 452 /* Jump operator P_D : results stored in pcis->vec1_B */ 453 454 ierr = VecScatterBegin(pcis->N_to_B,pcis->vec1_N,pcis->vec1_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 455 ierr = VecScatterEnd(pcis->N_to_B,pcis->vec1_N,pcis->vec1_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 456 /* Action of B_delta */ 457 ierr = MatMult(fetidpmat_ctx->B_delta,pcis->vec1_B,fetidpmat_ctx->lambda_local);CHKERRQ(ierr); 458 ierr = VecSet(fetidp_global,0.0);CHKERRQ(ierr); 459 ierr = VecScatterBegin(fetidpmat_ctx->l2g_lambda,fetidpmat_ctx->lambda_local,fetidp_global,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 460 ierr = VecScatterEnd(fetidpmat_ctx->l2g_lambda,fetidpmat_ctx->lambda_local,fetidp_global,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 461 /* Action of B_Ddelta^T */ 462 ierr = VecScatterBegin(fetidpmat_ctx->l2g_lambda,fetidp_global,fetidpmat_ctx->lambda_local,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 463 ierr = VecScatterEnd(fetidpmat_ctx->l2g_lambda,fetidp_global,fetidpmat_ctx->lambda_local,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 464 ierr = MatMultTranspose(fetidpmat_ctx->B_Ddelta,fetidpmat_ctx->lambda_local,pcis->vec1_B);CHKERRQ(ierr); 465 /* scaling */ 466 ierr = PCBDDCScalingExtension(fetidpmat_ctx->pc,pcis->vec1_B,pcis->vec1_global);CHKERRQ(ierr); 467 ierr = VecNorm(pcis->vec1_global,NORM_INFINITY,&val);CHKERRQ(ierr); 468 ierr = PetscViewerASCIIPrintf(viewer,"E: CHECK infty norm of R^T_D P_D: % 1.14e\n",val);CHKERRQ(ierr); 469 470 if (!fetidp->fully_redundant) { 471 /******************************************************************/ 472 /* TEST F: It should holds B_{delta}B^T_{D,delta}=I */ 473 /* Corollary thm 14 Mandel Tezaur and Dohrmann 2005 */ 474 /******************************************************************/ 475 ierr = VecDuplicate(fetidp_global,&test_vec);CHKERRQ(ierr); 476 ierr = VecSetRandom(fetidp_global,NULL);CHKERRQ(ierr); 477 if (fetidpmat_ctx->l2g_p) { 478 ierr = VecSet(fetidpmat_ctx->vP,0.);CHKERRQ(ierr); 479 ierr = VecScatterBegin(fetidpmat_ctx->l2g_p,fetidpmat_ctx->vP,fetidp_global,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 480 ierr = VecScatterEnd(fetidpmat_ctx->l2g_p,fetidpmat_ctx->vP,fetidp_global,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 481 } 482 /* Action of B_Ddelta^T */ 483 ierr = VecScatterBegin(fetidpmat_ctx->l2g_lambda,fetidp_global,fetidpmat_ctx->lambda_local,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 484 ierr = VecScatterEnd(fetidpmat_ctx->l2g_lambda,fetidp_global,fetidpmat_ctx->lambda_local,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 485 ierr = MatMultTranspose(fetidpmat_ctx->B_Ddelta,fetidpmat_ctx->lambda_local,pcis->vec1_B);CHKERRQ(ierr); 486 /* Action of B_delta */ 487 ierr = MatMult(fetidpmat_ctx->B_delta,pcis->vec1_B,fetidpmat_ctx->lambda_local);CHKERRQ(ierr); 488 ierr = VecSet(test_vec,0.0);CHKERRQ(ierr); 489 ierr = VecScatterBegin(fetidpmat_ctx->l2g_lambda,fetidpmat_ctx->lambda_local,test_vec,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 490 ierr = VecScatterEnd(fetidpmat_ctx->l2g_lambda,fetidpmat_ctx->lambda_local,test_vec,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 491 ierr = VecAXPY(fetidp_global,-1.,test_vec);CHKERRQ(ierr); 492 ierr = VecNorm(fetidp_global,NORM_INFINITY,&val);CHKERRQ(ierr); 493 ierr = PetscViewerASCIIPrintf(viewer,"E: CHECK infty norm of P^T_D - I: % 1.14e\n",val);CHKERRQ(ierr); 494 ierr = VecDestroy(&test_vec);CHKERRQ(ierr); 495 } 496 ierr = PetscViewerASCIIPrintf(viewer,"-------------------------------------\n");CHKERRQ(ierr); 497 ierr = PetscViewerFlush(viewer);CHKERRQ(ierr); 498 ierr = VecDestroy(&test_vec_p);CHKERRQ(ierr); 499 ierr = ISDestroy(&dirdofs);CHKERRQ(ierr); 500 ierr = VecDestroy(&fetidp_global);CHKERRQ(ierr); 501 ierr = ISRestoreIndices(isvert,&vertex_indices);CHKERRQ(ierr); 502 ierr = PCBDDCGraphRestoreCandidatesIS(pcbddc->mat_graph,NULL,NULL,NULL,NULL,&isvert);CHKERRQ(ierr); 503 PetscFunctionReturn(0); 504 } 505 506 #undef __FUNCT__ 507 #define __FUNCT__ "KSPFETIDPSetUpOperators" 508 static PetscErrorCode KSPFETIDPSetUpOperators(KSP ksp) 509 { 510 KSP_FETIDP *fetidp = (KSP_FETIDP*)ksp->data; 511 PC_BDDC *pcbddc = (PC_BDDC*)fetidp->innerbddc->data; 512 Mat A,Ap; 513 PetscInt fid = -1; 514 PetscBool ismatis,pisz; 515 PetscBool flip; /* Usually, Stokes is written 516 | A B'| | v | = | f | 517 | B 0 | | p | = | g | 518 If saddlepoint_flip is true, the code assumes it is written as 519 | A B'| | v | = | f | 520 |-B 0 | | p | = |-g | 521 */ 522 PetscObjectState matstate, matnnzstate; 523 PetscErrorCode ierr; 524 525 PetscFunctionBegin; 526 pisz = PETSC_FALSE; 527 flip = PETSC_FALSE; 528 ierr = PetscOptionsBegin(PetscObjectComm((PetscObject)ksp),((PetscObject)ksp)->prefix,"FETI-DP options","PC");CHKERRQ(ierr); 529 ierr = PetscOptionsInt("-ksp_fetidp_pressure_field","Field id for pressures for saddle-point problems",NULL,fid,&fid,NULL);CHKERRQ(ierr); 530 ierr = PetscOptionsBool("-ksp_fetidp_pressure_iszero","Zero pressure block",NULL,pisz,&pisz,NULL);CHKERRQ(ierr); 531 ierr = PetscOptionsBool("-ksp_fetidp_saddlepoint_flip","Flip the sign of the pressure-velocity (lower-left) block",NULL,flip,&flip,NULL);CHKERRQ(ierr); 532 ierr = PetscOptionsEnd();CHKERRQ(ierr); 533 534 fetidp->saddlepoint = (fid >= 0 ? PETSC_TRUE : fetidp->saddlepoint); 535 536 ierr = KSPGetOperators(ksp,&A,&Ap);CHKERRQ(ierr); 537 ierr = PetscObjectTypeCompare((PetscObject)Ap,MATIS,&ismatis);CHKERRQ(ierr); 538 if (!ismatis) SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_USER,"Pmat should be of type MATIS"); 539 540 /* Quiet return if the matrix states are unchanged. 541 Needed only for the saddle point case since it uses MatZeroRows 542 on a matrix that may not have changed */ 543 ierr = PetscObjectStateGet((PetscObject)Ap,&matstate);CHKERRQ(ierr); 544 ierr = MatGetNonzeroState(Ap,&matnnzstate);CHKERRQ(ierr); 545 if (matstate == fetidp->matstate && matnnzstate == fetidp->matnnzstate) PetscFunctionReturn(0); 546 fetidp->matstate = matstate; 547 fetidp->matnnzstate = matnnzstate; 548 fetidp->statechanged = fetidp->saddlepoint; 549 550 /* see if MATIS has same fields attached */ 551 if (!pcbddc->n_ISForDofsLocal && !pcbddc->n_ISForDofs) { 552 PetscContainer c; 553 554 ierr = PetscObjectQuery((PetscObject)Ap,"_convert_nest_lfields",(PetscObject*)&c);CHKERRQ(ierr); 555 if (c) { 556 MatISLocalFields lf; 557 ierr = PetscContainerGetPointer(c,(void**)&lf);CHKERRQ(ierr); 558 ierr = PCBDDCSetDofsSplittingLocal(fetidp->innerbddc,lf->nr,lf->rf);CHKERRQ(ierr); 559 } 560 } 561 562 if (!fetidp->saddlepoint) { 563 ierr = PCSetOperators(fetidp->innerbddc,A,Ap);CHKERRQ(ierr); 564 } else { 565 KSP kP; 566 Mat nA,lA; 567 Mat PAmat,PPmat; 568 IS pP; 569 PetscInt totP; 570 571 ierr = MatISGetLocalMat(Ap,&lA);CHKERRQ(ierr); 572 ierr = PetscObjectCompose((PetscObject)fetidp->innerbddc,"__KSPFETIDP_lA",(PetscObject)lA);CHKERRQ(ierr); 573 574 pP = fetidp->pP; 575 if (!pP) { /* first time, need to compute interface pressure dofs */ 576 PC_IS *pcis = (PC_IS*)fetidp->innerbddc->data; 577 Mat_IS *matis = (Mat_IS*)(Ap->data); 578 ISLocalToGlobalMapping l2g; 579 IS lP,II,pII,lPall,Pall,is1,is2; 580 const PetscInt *idxs; 581 PetscInt nl,ni,*widxs; 582 PetscInt i,j,n_neigh,*neigh,*n_shared,**shared,*count; 583 PetscInt rst,ren,n; 584 PetscBool ploc; 585 586 ierr = MatGetLocalSize(Ap,&nl,NULL);CHKERRQ(ierr); 587 ierr = MatGetOwnershipRange(Ap,&rst,&ren);CHKERRQ(ierr); 588 ierr = MatGetLocalSize(lA,&n,NULL);CHKERRQ(ierr); 589 590 if (!pcis->is_I_local) { /* need to compute interior dofs */ 591 ierr = PetscCalloc1(n,&count);CHKERRQ(ierr); 592 ierr = MatGetLocalToGlobalMapping(Ap,&l2g,NULL);CHKERRQ(ierr); 593 ierr = ISLocalToGlobalMappingGetInfo(l2g,&n_neigh,&neigh,&n_shared,&shared);CHKERRQ(ierr); 594 for (i=1;i<n_neigh;i++) 595 for (j=0;j<n_shared[i];j++) 596 count[shared[i][j]] += 1; 597 for (i=0,j=0;i<n;i++) if (!count[i]) count[j++] = i; 598 ierr = ISLocalToGlobalMappingRestoreInfo(l2g,&n_neigh,&neigh,&n_shared,&shared);CHKERRQ(ierr); 599 ierr = ISCreateGeneral(PETSC_COMM_SELF,j,count,PETSC_OWN_POINTER,&II);CHKERRQ(ierr); 600 } else { 601 ierr = PetscObjectReference((PetscObject)pcis->is_I_local);CHKERRQ(ierr); 602 II = pcis->is_I_local; 603 } 604 605 /* interior dofs in layout */ 606 ierr = MatISSetUpSF(Ap);CHKERRQ(ierr); 607 ierr = PetscMemzero(matis->sf_leafdata,n*sizeof(PetscInt));CHKERRQ(ierr); 608 ierr = PetscMemzero(matis->sf_rootdata,nl*sizeof(PetscInt));CHKERRQ(ierr); 609 ierr = ISGetLocalSize(II,&ni);CHKERRQ(ierr); 610 ierr = ISGetIndices(II,&idxs);CHKERRQ(ierr); 611 for (i=0;i<ni;i++) matis->sf_leafdata[idxs[i]] = 1; 612 ierr = ISRestoreIndices(II,&idxs);CHKERRQ(ierr); 613 ierr = PetscSFReduceBegin(matis->sf,MPIU_INT,matis->sf_leafdata,matis->sf_rootdata,MPIU_REPLACE);CHKERRQ(ierr); 614 ierr = PetscSFReduceEnd(matis->sf,MPIU_INT,matis->sf_leafdata,matis->sf_rootdata,MPIU_REPLACE);CHKERRQ(ierr); 615 ierr = PetscMalloc1(PetscMax(nl,n),&widxs);CHKERRQ(ierr); 616 for (i=0,ni=0;i<nl;i++) if (matis->sf_rootdata[i]) widxs[ni++] = i+rst; 617 ierr = ISCreateGeneral(PetscObjectComm((PetscObject)ksp),ni,widxs,PETSC_COPY_VALUES,&pII);CHKERRQ(ierr); 618 619 /* pressure space at the interface */ 620 Pall = NULL; 621 lPall = NULL; 622 ploc = PETSC_FALSE; 623 if (fid >= 0) { 624 if (pcbddc->n_ISForDofsLocal) { 625 PetscInt np; 626 627 if (fid >= pcbddc->n_ISForDofsLocal) SETERRQ2(PetscObjectComm((PetscObject)ksp),PETSC_ERR_USER,"Invalid field id for pressure %D, max %D",fid,pcbddc->n_ISForDofsLocal); 628 /* need a sequential IS */ 629 ierr = ISGetLocalSize(pcbddc->ISForDofsLocal[fid],&np);CHKERRQ(ierr); 630 ierr = ISGetIndices(pcbddc->ISForDofsLocal[fid],&idxs);CHKERRQ(ierr); 631 ierr = ISCreateGeneral(PETSC_COMM_SELF,np,idxs,PETSC_COPY_VALUES,&lPall);CHKERRQ(ierr); 632 ierr = ISRestoreIndices(pcbddc->ISForDofsLocal[fid],&idxs);CHKERRQ(ierr); 633 ploc = PETSC_TRUE; 634 } else if (pcbddc->n_ISForDofs) { 635 if (fid >= pcbddc->n_ISForDofs) SETERRQ2(PetscObjectComm((PetscObject)ksp),PETSC_ERR_USER,"Invalid field id for pressure %D, max %D",fid,pcbddc->n_ISForDofs); 636 ierr = PetscObjectReference((PetscObject)pcbddc->ISForDofs[fid]);CHKERRQ(ierr); 637 Pall = pcbddc->ISForDofs[fid]; 638 } else SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_USER,"Missing fields! Use PCBDDCSetDofsSplitting/Local"); 639 } else { /* fallback to zero pressure block */ 640 IS list[2]; 641 642 ierr = MatFindZeroDiagonals(Ap,&list[1]);CHKERRQ(ierr); 643 ierr = ISComplement(list[1],rst,ren,&list[0]);CHKERRQ(ierr); 644 ierr = PCBDDCSetDofsSplitting(fetidp->innerbddc,2,list);CHKERRQ(ierr); 645 ierr = ISDestroy(&list[0]);CHKERRQ(ierr); 646 Pall = list[1]; 647 } 648 if (ploc) { 649 ierr = ISDifference(lPall,II,&lP);CHKERRQ(ierr); 650 ierr = PetscObjectCompose((PetscObject)fetidp->innerbddc,"__KSPFETIDP_lP",(PetscObject)lP);CHKERRQ(ierr); 651 } else { 652 ierr = ISDifference(Pall,pII,&pP);CHKERRQ(ierr); 653 ierr = PetscObjectCompose((PetscObject)fetidp->innerbddc,"__KSPFETIDP_pP",(PetscObject)pP);CHKERRQ(ierr); 654 /* need all local pressure dofs */ 655 ierr = PetscMemzero(matis->sf_leafdata,n*sizeof(PetscInt));CHKERRQ(ierr); 656 ierr = PetscMemzero(matis->sf_rootdata,nl*sizeof(PetscInt));CHKERRQ(ierr); 657 ierr = ISGetLocalSize(Pall,&ni);CHKERRQ(ierr); 658 ierr = ISGetIndices(Pall,&idxs);CHKERRQ(ierr); 659 for (i=0;i<ni;i++) matis->sf_rootdata[idxs[i]-rst] = 1; 660 ierr = ISRestoreIndices(Pall,&idxs);CHKERRQ(ierr); 661 ierr = PetscSFBcastBegin(matis->sf,MPIU_INT,matis->sf_rootdata,matis->sf_leafdata);CHKERRQ(ierr); 662 ierr = PetscSFBcastEnd(matis->sf,MPIU_INT,matis->sf_rootdata,matis->sf_leafdata);CHKERRQ(ierr); 663 for (i=0,ni=0;i<n;i++) if (matis->sf_leafdata[i]) widxs[ni++] = i; 664 ierr = ISCreateGeneral(PETSC_COMM_SELF,ni,widxs,PETSC_COPY_VALUES,&lPall);CHKERRQ(ierr); 665 } 666 667 if (flip) { 668 PetscInt npl; 669 if (!Pall) { 670 ierr = PetscMemzero(matis->sf_leafdata,n*sizeof(PetscInt));CHKERRQ(ierr); 671 ierr = PetscMemzero(matis->sf_rootdata,nl*sizeof(PetscInt));CHKERRQ(ierr); 672 ierr = ISGetLocalSize(lPall,&ni);CHKERRQ(ierr); 673 ierr = ISGetIndices(lPall,&idxs);CHKERRQ(ierr); 674 for (i=0;i<ni;i++) matis->sf_leafdata[idxs[i]] = 1; 675 ierr = ISRestoreIndices(lPall,&idxs);CHKERRQ(ierr); 676 ierr = PetscSFReduceBegin(matis->sf,MPIU_INT,matis->sf_leafdata,matis->sf_rootdata,MPIU_REPLACE);CHKERRQ(ierr); 677 ierr = PetscSFReduceEnd(matis->sf,MPIU_INT,matis->sf_leafdata,matis->sf_rootdata,MPIU_REPLACE);CHKERRQ(ierr); 678 for (i=0,ni=0;i<nl;i++) if (matis->sf_rootdata[i]) widxs[ni++] = i+rst; 679 ierr = ISCreateGeneral(PetscObjectComm((PetscObject)ksp),ni,widxs,PETSC_COPY_VALUES,&Pall);CHKERRQ(ierr); 680 } 681 ierr = ISGetLocalSize(Pall,&npl);CHKERRQ(ierr); 682 ierr = ISGetIndices(Pall,&idxs);CHKERRQ(ierr); 683 ierr = MatCreateVecs(Ap,NULL,&fetidp->rhs_flip);CHKERRQ(ierr); 684 ierr = VecSet(fetidp->rhs_flip,1.);CHKERRQ(ierr); 685 ierr = VecSetOption(fetidp->rhs_flip,VEC_IGNORE_OFF_PROC_ENTRIES,PETSC_TRUE);CHKERRQ(ierr); 686 for (i=0;i<npl;i++) { 687 ierr = VecSetValue(fetidp->rhs_flip,idxs[i],-1.,INSERT_VALUES);CHKERRQ(ierr); 688 } 689 ierr = VecAssemblyBegin(fetidp->rhs_flip);CHKERRQ(ierr); 690 ierr = VecAssemblyEnd(fetidp->rhs_flip);CHKERRQ(ierr); 691 ierr = PetscObjectCompose((PetscObject)fetidp->innerbddc,"__KSPFETIDP_flip",(PetscObject)fetidp->rhs_flip);CHKERRQ(ierr); 692 ierr = ISRestoreIndices(Pall,&idxs);CHKERRQ(ierr); 693 } 694 ierr = ISDestroy(&Pall);CHKERRQ(ierr); 695 ierr = ISDestroy(&pII);CHKERRQ(ierr); 696 697 /* local interface pressures in subdomain-wise and global ordering */ 698 ierr = PetscMemzero(matis->sf_leafdata,n*sizeof(PetscInt));CHKERRQ(ierr); 699 ierr = PetscMemzero(matis->sf_rootdata,nl*sizeof(PetscInt));CHKERRQ(ierr); 700 if (pP) { 701 ierr = ISGetLocalSize(pP,&ni);CHKERRQ(ierr); 702 ierr = ISGetIndices(pP,&idxs);CHKERRQ(ierr); 703 for (i=0;i<ni;i++) matis->sf_rootdata[idxs[i]-rst] = 1; 704 ierr = ISRestoreIndices(pP,&idxs);CHKERRQ(ierr); 705 ierr = PetscSFBcastBegin(matis->sf,MPIU_INT,matis->sf_rootdata,matis->sf_leafdata);CHKERRQ(ierr); 706 ierr = PetscSFBcastEnd(matis->sf,MPIU_INT,matis->sf_rootdata,matis->sf_leafdata);CHKERRQ(ierr); 707 for (i=0,ni=0;i<n;i++) if (matis->sf_leafdata[i]) widxs[ni++] = i; 708 ierr = ISLocalToGlobalMappingApply(l2g,ni,widxs,widxs+ni);CHKERRQ(ierr); 709 ierr = ISCreateGeneral(PETSC_COMM_SELF,ni,widxs,PETSC_COPY_VALUES,&lP);CHKERRQ(ierr); 710 ierr = PetscObjectCompose((PetscObject)fetidp->innerbddc,"__KSPFETIDP_lP",(PetscObject)lP);CHKERRQ(ierr); 711 ierr = ISCreateGeneral(PetscObjectComm((PetscObject)ksp),ni,widxs+ni,PETSC_COPY_VALUES,&is1);CHKERRQ(ierr); 712 ierr = PetscObjectCompose((PetscObject)fetidp->innerbddc,"__KSPFETIDP_gP",(PetscObject)is1);CHKERRQ(ierr); 713 ierr = ISDestroy(&is1);CHKERRQ(ierr); 714 } else { 715 ierr = ISGetLocalSize(lP,&ni);CHKERRQ(ierr); 716 ierr = ISGetIndices(lP,&idxs);CHKERRQ(ierr); 717 for (i=0;i<ni;i++) 718 if (idxs[i] >=0 && idxs[i] < n) 719 matis->sf_leafdata[idxs[i]] = 1; 720 ierr = ISRestoreIndices(lP,&idxs);CHKERRQ(ierr); 721 ierr = PetscSFReduceBegin(matis->sf,MPIU_INT,matis->sf_leafdata,matis->sf_rootdata,MPIU_REPLACE);CHKERRQ(ierr); 722 ierr = ISLocalToGlobalMappingApply(l2g,ni,idxs,widxs);CHKERRQ(ierr); 723 ierr = ISCreateGeneral(PetscObjectComm((PetscObject)ksp),ni,widxs,PETSC_COPY_VALUES,&is1);CHKERRQ(ierr); 724 ierr = PetscObjectCompose((PetscObject)fetidp->innerbddc,"__KSPFETIDP_gP",(PetscObject)is1);CHKERRQ(ierr); 725 ierr = ISDestroy(&is1);CHKERRQ(ierr); 726 ierr = PetscSFReduceEnd(matis->sf,MPIU_INT,matis->sf_leafdata,matis->sf_rootdata,MPIU_REPLACE);CHKERRQ(ierr); 727 for (i=0,ni=0;i<nl;i++) if (matis->sf_rootdata[i]) widxs[ni++] = i+rst; 728 ierr = ISCreateGeneral(PetscObjectComm((PetscObject)ksp),ni,widxs,PETSC_COPY_VALUES,&pP);CHKERRQ(ierr); 729 ierr = PetscObjectCompose((PetscObject)fetidp->innerbddc,"__KSPFETIDP_pP",(PetscObject)pP);CHKERRQ(ierr); 730 } 731 ierr = PetscFree(widxs);CHKERRQ(ierr); 732 733 /* We may want to use a discrete harmonic solver instead 734 of a Stokes harmonic for the Dirichler preconditioner 735 Need to extract the interior pressure dofs in interior dofs ordering */ 736 ierr = ISDifference(lPall,lP,&is1);CHKERRQ(ierr); 737 ierr = ISDifference(II,is1,&is2);CHKERRQ(ierr); 738 ierr = ISDestroy(&is1);CHKERRQ(ierr); 739 ierr = ISLocalToGlobalMappingCreateIS(II,&l2g);CHKERRQ(ierr); 740 ierr = ISGlobalToLocalMappingApplyIS(l2g,IS_GTOLM_DROP,is2,&is1);CHKERRQ(ierr); 741 ierr = ISGetLocalSize(is1,&i);CHKERRQ(ierr); 742 ierr = ISGetLocalSize(is2,&j);CHKERRQ(ierr); 743 if (i != j) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Inconsistent local sizes %D and %D for iP",i,j); 744 ierr = PetscObjectCompose((PetscObject)fetidp->innerbddc,"__KSPFETIDP_iP",(PetscObject)is1);CHKERRQ(ierr); 745 ierr = ISLocalToGlobalMappingDestroy(&l2g);CHKERRQ(ierr); 746 ierr = ISDestroy(&is1);CHKERRQ(ierr); 747 ierr = ISDestroy(&is2);CHKERRQ(ierr); 748 ierr = ISDestroy(&lPall);CHKERRQ(ierr); 749 ierr = ISDestroy(&II);CHKERRQ(ierr); 750 751 /* exclude interface pressures from the inner BDDC */ 752 if (pcbddc->DirichletBoundariesLocal) { 753 IS list[2],plP,isout; 754 PetscInt np; 755 756 /* need a parallel IS */ 757 ierr = ISGetLocalSize(lP,&np);CHKERRQ(ierr); 758 ierr = ISGetIndices(lP,&idxs);CHKERRQ(ierr); 759 ierr = ISCreateGeneral(PetscObjectComm((PetscObject)ksp),np,idxs,PETSC_USE_POINTER,&plP);CHKERRQ(ierr); 760 list[0] = plP; 761 list[1] = pcbddc->DirichletBoundariesLocal; 762 ierr = ISConcatenate(PetscObjectComm((PetscObject)ksp),2,list,&isout);CHKERRQ(ierr); 763 ierr = ISDestroy(&plP);CHKERRQ(ierr); 764 ierr = ISRestoreIndices(lP,&idxs);CHKERRQ(ierr); 765 ierr = PCBDDCSetDirichletBoundariesLocal(fetidp->innerbddc,isout);CHKERRQ(ierr); 766 ierr = ISDestroy(&isout);CHKERRQ(ierr); 767 } else if (pcbddc->DirichletBoundaries) { 768 IS list[2],isout; 769 770 list[0] = pP; 771 list[1] = pcbddc->DirichletBoundaries; 772 ierr = ISConcatenate(PetscObjectComm((PetscObject)ksp),2,list,&isout);CHKERRQ(ierr); 773 ierr = PCBDDCSetDirichletBoundaries(fetidp->innerbddc,isout);CHKERRQ(ierr); 774 ierr = ISDestroy(&isout);CHKERRQ(ierr); 775 } else { 776 IS plP; 777 PetscInt np; 778 779 /* need a parallel IS */ 780 ierr = ISGetLocalSize(lP,&np);CHKERRQ(ierr); 781 ierr = ISGetIndices(lP,&idxs);CHKERRQ(ierr); 782 ierr = ISCreateGeneral(PetscObjectComm((PetscObject)ksp),np,idxs,PETSC_COPY_VALUES,&plP);CHKERRQ(ierr); 783 ierr = PCBDDCSetDirichletBoundariesLocal(fetidp->innerbddc,plP);CHKERRQ(ierr); 784 ierr = ISDestroy(&plP);CHKERRQ(ierr); 785 ierr = ISRestoreIndices(lP,&idxs);CHKERRQ(ierr); 786 } 787 ierr = ISDestroy(&lP);CHKERRQ(ierr); 788 fetidp->pP = pP; 789 } 790 791 /* total number of pressure interface dofs */ 792 ierr = ISGetSize(fetidp->pP,&totP);CHKERRQ(ierr); 793 794 /* Set operator for inner BDDC */ 795 if (totP || fetidp->rhs_flip) { 796 ierr = MatDuplicate(Ap,MAT_COPY_VALUES,&nA);CHKERRQ(ierr); 797 } else { 798 ierr = PetscObjectReference((PetscObject)Ap);CHKERRQ(ierr); 799 nA = Ap; 800 } 801 if (fetidp->rhs_flip) { 802 Mat lA2; 803 804 ierr = MatDiagonalScale(nA,fetidp->rhs_flip,NULL);CHKERRQ(ierr); 805 if (totP) { 806 ierr = MatISGetLocalMat(nA,&lA);CHKERRQ(ierr); 807 ierr = MatDuplicate(lA,MAT_COPY_VALUES,&lA2);CHKERRQ(ierr); 808 ierr = PetscObjectCompose((PetscObject)fetidp->innerbddc,"__KSPFETIDP_lA",(PetscObject)lA2);CHKERRQ(ierr); 809 } 810 ierr = MatDestroy(&lA2);CHKERRQ(ierr); 811 } 812 if (totP) { 813 ierr = MatSetOption(nA,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_FALSE);CHKERRQ(ierr); 814 ierr = MatZeroRowsColumnsIS(nA,fetidp->pP,1.,NULL,NULL);CHKERRQ(ierr); 815 } else { 816 ierr = PetscObjectCompose((PetscObject)fetidp->innerbddc,"__KSPFETIDP_lA",NULL);CHKERRQ(ierr); 817 } 818 ierr = PCSetOperators(fetidp->innerbddc,nA,nA);CHKERRQ(ierr); 819 ierr = MatDestroy(&nA);CHKERRQ(ierr); 820 821 /* non-zero rhs on interior dofs when applying the preconditioner */ 822 if (totP) pcbddc->switch_static = PETSC_TRUE; 823 824 /* if there are no interface pressures, set inner bddc flag for benign saddle point */ 825 if (!totP) pcbddc->benign_saddle_point = PETSC_TRUE; 826 827 /* Operators for pressure preconditioner */ 828 if (totP) { 829 ierr = PetscObjectQuery((PetscObject)fetidp->innerbddc,"__KSPFETIDP_PAmat",(PetscObject*)&PAmat);CHKERRQ(ierr); 830 ierr = PetscObjectQuery((PetscObject)fetidp->innerbddc,"__KSPFETIDP_PPmat",(PetscObject*)&PPmat);CHKERRQ(ierr); 831 if (PAmat) { 832 ierr = PetscObjectReference((PetscObject)PAmat);CHKERRQ(ierr); 833 } 834 if (PPmat) { 835 ierr = PetscObjectReference((PetscObject)PPmat);CHKERRQ(ierr); 836 } 837 838 /* Extract pressure block */ 839 if (!pisz) { 840 Mat C; 841 842 ierr = MatGetSubMatrix(Ap,fetidp->pP,fetidp->pP,MAT_INITIAL_MATRIX,&C);CHKERRQ(ierr); 843 ierr = MatScale(C,-1.);CHKERRQ(ierr); 844 ierr = PetscObjectCompose((PetscObject)fetidp->innerbddc,"__KSPFETIDP_C",(PetscObject)C);CHKERRQ(ierr); 845 /* default operators for the interface pressure solver */ 846 if (!PAmat) { 847 ierr = PetscObjectReference((PetscObject)C);CHKERRQ(ierr); 848 PAmat = C; 849 ierr = PetscObjectCompose((PetscObject)fetidp->innerbddc,"__KSPFETIDP_PAmat",(PetscObject)PAmat);CHKERRQ(ierr); 850 } 851 if (!PPmat) { 852 ierr = PetscObjectReference((PetscObject)C);CHKERRQ(ierr); 853 PPmat = C; 854 ierr = PetscObjectCompose((PetscObject)fetidp->innerbddc,"__KSPFETIDP_PPmat",(PetscObject)PPmat);CHKERRQ(ierr); 855 } 856 ierr = MatDestroy(&C);CHKERRQ(ierr); 857 } 858 859 if (!PAmat) { /* Just use the identity */ 860 PetscInt nl; 861 862 ierr = ISGetLocalSize(fetidp->pP,&nl);CHKERRQ(ierr); 863 ierr = MatCreate(PetscObjectComm((PetscObject)ksp),&PAmat);CHKERRQ(ierr); 864 ierr = MatSetSizes(PAmat,nl,nl,totP,totP);CHKERRQ(ierr); 865 ierr = MatSetType(PAmat,MATAIJ);CHKERRQ(ierr); 866 ierr = MatMPIAIJSetPreallocation(PAmat,1,NULL,0,NULL);CHKERRQ(ierr); 867 ierr = MatSeqAIJSetPreallocation(PAmat,1,NULL);CHKERRQ(ierr); 868 ierr = MatAssemblyBegin(PAmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 869 ierr = MatAssemblyEnd(PAmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 870 ierr = MatShift(PAmat,1.);CHKERRQ(ierr); 871 ierr = PetscObjectCompose((PetscObject)fetidp->innerbddc,"__KSPFETIDP_PAmat",(PetscObject)PAmat);CHKERRQ(ierr); 872 } else { /* check size of pressure operator and restrict if needed */ 873 PetscInt AM,PAM,PAN,pam,pan,am,an,pl; 874 875 ierr = MatGetSize(Ap,&AM,NULL);CHKERRQ(ierr); 876 ierr = MatGetSize(PAmat,&PAM,&PAN);CHKERRQ(ierr); 877 ierr = MatGetLocalSize(PAmat,&pam,&pan);CHKERRQ(ierr); 878 ierr = MatGetLocalSize(Ap,&am,&an);CHKERRQ(ierr); 879 ierr = ISGetLocalSize(fetidp->pP,&pl);CHKERRQ(ierr); 880 if (PAM != PAN) SETERRQ2(PetscObjectComm((PetscObject)ksp),PETSC_ERR_USER,"Pressure matrix must be square, unsupported %D x %D",PAM,PAN); 881 if (pam != pan) SETERRQ2(PetscObjectComm((PetscObject)ksp),PETSC_ERR_USER,"Local sizes of pressure matrix must be equal, unsupported %D x %D",pam,pan); 882 if (pam != am && pam != pl) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_USER,"Invalid number of local rows %D for pressure matrix! Supported are %D or %D",pam,am,pl); 883 if (pan != an && pan != pl) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_USER,"Invalid number of local columns %D for pressure matrix! Supported are %D or %D",pan,an,pl); 884 if (PAM == AM) { /* monolithic ordering, restrict to interface pressure */ 885 Mat C; 886 887 ierr = MatGetSubMatrix(PAmat,fetidp->pP,fetidp->pP,MAT_INITIAL_MATRIX,&C);CHKERRQ(ierr); 888 ierr = MatDestroy(&PAmat);CHKERRQ(ierr); 889 PAmat = C; 890 ierr = PetscObjectCompose((PetscObject)fetidp->innerbddc,"__KSPFETIDP_PAmat",(PetscObject)PAmat);CHKERRQ(ierr); 891 } 892 } 893 if (PPmat) { 894 PetscInt AM,PAM,PAN,pam,pan,am,an,pl; 895 896 ierr = MatGetSize(Ap,&AM,NULL);CHKERRQ(ierr); 897 ierr = MatGetSize(PPmat,&PAM,&PAN);CHKERRQ(ierr); 898 ierr = MatGetLocalSize(PPmat,&pam,&pan);CHKERRQ(ierr); 899 ierr = MatGetLocalSize(Ap,&am,&an);CHKERRQ(ierr); 900 ierr = ISGetLocalSize(fetidp->pP,&pl);CHKERRQ(ierr); 901 if (PAM != PAN) SETERRQ2(PetscObjectComm((PetscObject)ksp),PETSC_ERR_USER,"Pressure matrix must be square, unsupported %D x %D",PAM,PAN); 902 if (pam != pan) SETERRQ2(PetscObjectComm((PetscObject)ksp),PETSC_ERR_USER,"Local sizes of pressure matrix must be equal, unsupported %D x %D",pam,pan); 903 if (pam != am && pam != pl) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_USER,"Invalid number of local rows %D for pressure matrix! Supported are %D or %D",pam,am,pl); 904 if (pan != an && pan != pl) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_USER,"Invalid number of local columns %D for pressure matrix! Supported are %D or %D",pan,an,pl); 905 if (PAM == AM) { /* monolithic ordering, restrict to interface pressure */ 906 Mat C; 907 908 ierr = MatGetSubMatrix(PPmat,fetidp->pP,fetidp->pP,MAT_INITIAL_MATRIX,&C);CHKERRQ(ierr); 909 ierr = MatDestroy(&PPmat);CHKERRQ(ierr); 910 PPmat = C; 911 ierr = PetscObjectCompose((PetscObject)fetidp->innerbddc,"__KSPFETIDP_PPmat",(PetscObject)PPmat);CHKERRQ(ierr); 912 } 913 } else { 914 ierr = PetscObjectReference((PetscObject)PAmat);CHKERRQ(ierr); 915 PPmat = PAmat; 916 } 917 918 /* create pressure solver */ 919 ierr = KSPCreate(PetscObjectComm((PetscObject)ksp),&kP);CHKERRQ(ierr); 920 ierr = KSPSetOperators(kP,PAmat,PPmat);CHKERRQ(ierr); 921 ierr = KSPSetOptionsPrefix(kP,((PetscObject)ksp)->prefix);CHKERRQ(ierr); 922 ierr = KSPAppendOptionsPrefix(kP,"fetidp_p_");CHKERRQ(ierr); 923 ierr = KSPSetFromOptions(kP);CHKERRQ(ierr); 924 ierr = PetscObjectCompose((PetscObject)fetidp->innerbddc,"__KSPFETIDP_PKSP",(PetscObject)kP);CHKERRQ(ierr); 925 ierr = MatDestroy(&PAmat);CHKERRQ(ierr); 926 ierr = MatDestroy(&PPmat);CHKERRQ(ierr); 927 ierr = KSPDestroy(&kP);CHKERRQ(ierr); 928 } else { /* totP == 0 */ 929 ierr = PetscObjectCompose((PetscObject)fetidp->innerbddc,"__KSPFETIDP_pP",NULL);CHKERRQ(ierr); 930 } 931 } 932 PetscFunctionReturn(0); 933 } 934 935 #undef __FUNCT__ 936 #define __FUNCT__ "KSPSetUp_FETIDP" 937 static PetscErrorCode KSPSetUp_FETIDP(KSP ksp) 938 { 939 KSP_FETIDP *fetidp = (KSP_FETIDP*)ksp->data; 940 PC_BDDC *pcbddc = (PC_BDDC*)fetidp->innerbddc->data; 941 PetscBool flg; 942 PetscErrorCode ierr; 943 944 PetscFunctionBegin; 945 ierr = KSPFETIDPSetUpOperators(ksp);CHKERRQ(ierr); 946 /* set up BDDC */ 947 ierr = PCSetUp(fetidp->innerbddc);CHKERRQ(ierr); 948 949 /* FETI-DP as it is implemented needs an exact coarse solver */ 950 if (pcbddc->coarse_ksp) { 951 ierr = KSPSetTolerances(pcbddc->coarse_ksp,PETSC_SMALL,PETSC_SMALL,PETSC_DEFAULT,1000);CHKERRQ(ierr); 952 ierr = KSPSetNormType(pcbddc->coarse_ksp,KSP_NORM_DEFAULT);CHKERRQ(ierr); 953 } 954 /* FETI-DP as it is implemented needs exact local Neumann solvers */ 955 ierr = KSPSetTolerances(pcbddc->ksp_R,PETSC_SMALL,PETSC_SMALL,PETSC_DEFAULT,1000);CHKERRQ(ierr); 956 ierr = KSPSetNormType(pcbddc->ksp_R,KSP_NORM_DEFAULT);CHKERRQ(ierr); 957 958 /* setup FETI-DP operators 959 If fetidp->statechanged is true, we need update the operators 960 that are needed in the saddle-point case. This should be replaced 961 by a better logic when the FETI-DP matrix and preconditioner will 962 have their own classes */ 963 if (pcbddc->new_primal_space || fetidp->statechanged) { 964 Mat F; /* the FETI-DP matrix */ 965 PC D; /* the FETI-DP preconditioner */ 966 ierr = KSPReset(fetidp->innerksp);CHKERRQ(ierr); 967 ierr = PCBDDCCreateFETIDPOperators(fetidp->innerbddc,fetidp->fully_redundant,((PetscObject)ksp)->prefix,&F,&D);CHKERRQ(ierr); 968 ierr = KSPSetOperators(fetidp->innerksp,F,F);CHKERRQ(ierr); 969 ierr = KSPSetTolerances(fetidp->innerksp,ksp->rtol,ksp->abstol,ksp->divtol,ksp->max_it);CHKERRQ(ierr); 970 ierr = KSPSetPC(fetidp->innerksp,D);CHKERRQ(ierr); 971 ierr = KSPSetFromOptions(fetidp->innerksp);CHKERRQ(ierr); 972 ierr = MatCreateVecs(F,&(fetidp->innerksp)->vec_rhs,&(fetidp->innerksp)->vec_sol);CHKERRQ(ierr); 973 ierr = MatDestroy(&F);CHKERRQ(ierr); 974 ierr = PCDestroy(&D);CHKERRQ(ierr); 975 if (fetidp->check) { 976 PetscViewer viewer; 977 978 if (!pcbddc->dbg_viewer) { 979 viewer = PETSC_VIEWER_STDOUT_(PetscObjectComm((PetscObject)ksp)); 980 } else { 981 viewer = pcbddc->dbg_viewer; 982 } 983 ierr = KSPFETIDPCheckOperators(ksp,viewer);CHKERRQ(ierr); 984 } 985 } 986 fetidp->statechanged = PETSC_FALSE; 987 pcbddc->new_primal_space = PETSC_FALSE; 988 989 /* propagate settings to the inner solve */ 990 ierr = KSPGetComputeSingularValues(ksp,&flg);CHKERRQ(ierr); 991 ierr = KSPSetComputeSingularValues(fetidp->innerksp,flg);CHKERRQ(ierr); 992 if (ksp->res_hist) { 993 ierr = KSPSetResidualHistory(fetidp->innerksp,ksp->res_hist,ksp->res_hist_max,ksp->res_hist_reset);CHKERRQ(ierr); 994 } 995 ierr = KSPSetUp(fetidp->innerksp);CHKERRQ(ierr); 996 PetscFunctionReturn(0); 997 } 998 999 #undef __FUNCT__ 1000 #define __FUNCT__ "KSPSolve_FETIDP" 1001 static PetscErrorCode KSPSolve_FETIDP(KSP ksp) 1002 { 1003 PetscErrorCode ierr; 1004 Mat F,A; 1005 MatNullSpace nsp; 1006 Vec X,B,Xl,Bl; 1007 KSP_FETIDP *fetidp = (KSP_FETIDP*)ksp->data; 1008 PC_BDDC *pcbddc = (PC_BDDC*)fetidp->innerbddc->data; 1009 1010 PetscFunctionBegin; 1011 ierr = KSPGetOperators(ksp,&A,NULL);CHKERRQ(ierr); 1012 ierr = KSPGetRhs(ksp,&B);CHKERRQ(ierr); 1013 ierr = KSPGetSolution(ksp,&X);CHKERRQ(ierr); 1014 ierr = KSPGetOperators(fetidp->innerksp,&F,NULL);CHKERRQ(ierr); 1015 ierr = KSPGetRhs(fetidp->innerksp,&Bl);CHKERRQ(ierr); 1016 ierr = KSPGetSolution(fetidp->innerksp,&Xl);CHKERRQ(ierr); 1017 ierr = PCBDDCMatFETIDPGetRHS(F,B,Bl);CHKERRQ(ierr); 1018 if (ksp->transpose_solve) { 1019 ierr = KSPSolveTranspose(fetidp->innerksp,Bl,Xl);CHKERRQ(ierr); 1020 } else { 1021 ierr = KSPSolve(fetidp->innerksp,Bl,Xl);CHKERRQ(ierr); 1022 } 1023 ierr = PCBDDCMatFETIDPGetSolution(F,Xl,X);CHKERRQ(ierr); 1024 ierr = MatGetNullSpace(A,&nsp);CHKERRQ(ierr); 1025 if (nsp) { 1026 ierr = MatNullSpaceRemove(nsp,X);CHKERRQ(ierr); 1027 } 1028 /* update ksp with stats from inner ksp */ 1029 ierr = KSPGetConvergedReason(fetidp->innerksp,&ksp->reason);CHKERRQ(ierr); 1030 ierr = KSPGetIterationNumber(fetidp->innerksp,&ksp->its);CHKERRQ(ierr); 1031 ksp->totalits += ksp->its; 1032 ierr = KSPGetResidualHistory(fetidp->innerksp,NULL,&ksp->res_hist_len);CHKERRQ(ierr); 1033 /* restore defaults for inner BDDC (Pre/PostSolve flags) */ 1034 pcbddc->temp_solution_used = PETSC_FALSE; 1035 pcbddc->rhs_change = PETSC_FALSE; 1036 pcbddc->exact_dirichlet_trick_app = PETSC_FALSE; 1037 PetscFunctionReturn(0); 1038 } 1039 1040 #undef __FUNCT__ 1041 #define __FUNCT__ "KSPReset_FETIDP" 1042 static PetscErrorCode KSPReset_FETIDP(KSP ksp) 1043 { 1044 KSP_FETIDP *fetidp = (KSP_FETIDP*)ksp->data; 1045 PetscErrorCode ierr; 1046 1047 PetscFunctionBegin; 1048 ierr = ISDestroy(&fetidp->pP);CHKERRQ(ierr); 1049 ierr = VecDestroy(&fetidp->rhs_flip);CHKERRQ(ierr); 1050 ierr = PCReset(fetidp->innerbddc);CHKERRQ(ierr); 1051 ierr = KSPReset(fetidp->innerksp);CHKERRQ(ierr); 1052 fetidp->saddlepoint = PETSC_FALSE; 1053 fetidp->matstate = -1; 1054 fetidp->matnnzstate = -1; 1055 fetidp->statechanged = PETSC_TRUE; 1056 PetscFunctionReturn(0); 1057 } 1058 1059 #undef __FUNCT__ 1060 #define __FUNCT__ "KSPDestroy_FETIDP" 1061 static PetscErrorCode KSPDestroy_FETIDP(KSP ksp) 1062 { 1063 KSP_FETIDP *fetidp = (KSP_FETIDP*)ksp->data; 1064 PetscErrorCode ierr; 1065 1066 PetscFunctionBegin; 1067 ierr = KSPReset_FETIDP(ksp);CHKERRQ(ierr); 1068 ierr = PCDestroy(&fetidp->innerbddc);CHKERRQ(ierr); 1069 ierr = KSPDestroy(&fetidp->innerksp);CHKERRQ(ierr); 1070 ierr = PetscFree(fetidp->monctx);CHKERRQ(ierr); 1071 ierr = PetscObjectComposeFunction((PetscObject)ksp,"KSPFETIDPSetInnerBDDC_C",NULL);CHKERRQ(ierr); 1072 ierr = PetscObjectComposeFunction((PetscObject)ksp,"KSPFETIDPGetInnerBDDC_C",NULL);CHKERRQ(ierr); 1073 ierr = PetscObjectComposeFunction((PetscObject)ksp,"KSPFETIDPGetInnerKSP_C",NULL);CHKERRQ(ierr); 1074 ierr = PetscObjectComposeFunction((PetscObject)ksp,"KSPFETIDPSetPressureOperators_C",NULL);CHKERRQ(ierr); 1075 ierr = PetscFree(ksp->data);CHKERRQ(ierr); 1076 PetscFunctionReturn(0); 1077 } 1078 1079 #undef __FUNCT__ 1080 #define __FUNCT__ "KSPView_FETIDP" 1081 static PetscErrorCode KSPView_FETIDP(KSP ksp,PetscViewer viewer) 1082 { 1083 KSP_FETIDP *fetidp = (KSP_FETIDP*)ksp->data; 1084 PetscErrorCode ierr; 1085 PetscBool iascii; 1086 1087 PetscFunctionBegin; 1088 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr); 1089 if (iascii) { 1090 ierr = PetscViewerASCIIPrintf(viewer," FETI-DP: fully redundant: %D\n",fetidp->fully_redundant);CHKERRQ(ierr); 1091 ierr = PetscViewerASCIIPrintf(viewer," FETI-DP: saddle point: %D\n",fetidp->saddlepoint);CHKERRQ(ierr); 1092 ierr = PetscViewerASCIIPrintf(viewer," FETI-DP: inner solver details\n");CHKERRQ(ierr); 1093 ierr = PetscViewerASCIIAddTab(viewer,2);CHKERRQ(ierr); 1094 } 1095 ierr = KSPView(fetidp->innerksp,viewer);CHKERRQ(ierr); 1096 if (iascii) { 1097 ierr = PetscViewerASCIISubtractTab(viewer,2);CHKERRQ(ierr); 1098 ierr = PetscViewerASCIIPrintf(viewer," FETI-DP: BDDC solver details\n");CHKERRQ(ierr); 1099 ierr = PetscViewerASCIIAddTab(viewer,2);CHKERRQ(ierr); 1100 } 1101 ierr = PCView(fetidp->innerbddc,viewer);CHKERRQ(ierr); 1102 if (iascii) { 1103 ierr = PetscViewerASCIISubtractTab(viewer,2);CHKERRQ(ierr); 1104 } 1105 PetscFunctionReturn(0); 1106 } 1107 1108 #undef __FUNCT__ 1109 #define __FUNCT__ "KSPSetFromOptions_FETIDP" 1110 static PetscErrorCode KSPSetFromOptions_FETIDP(PetscOptionItems *PetscOptionsObject,KSP ksp) 1111 { 1112 KSP_FETIDP *fetidp = (KSP_FETIDP*)ksp->data; 1113 PetscErrorCode ierr; 1114 1115 PetscFunctionBegin; 1116 /* set options prefixes for the inner objects, since the parent prefix will be valid at this point */ 1117 ierr = PetscObjectSetOptionsPrefix((PetscObject)fetidp->innerksp,((PetscObject)ksp)->prefix);CHKERRQ(ierr); 1118 ierr = PetscObjectAppendOptionsPrefix((PetscObject)fetidp->innerksp,"fetidp_");CHKERRQ(ierr); 1119 if (!fetidp->userbddc) { 1120 ierr = PetscObjectSetOptionsPrefix((PetscObject)fetidp->innerbddc,((PetscObject)ksp)->prefix);CHKERRQ(ierr); 1121 ierr = PetscObjectAppendOptionsPrefix((PetscObject)fetidp->innerbddc,"fetidp_bddc_");CHKERRQ(ierr); 1122 } 1123 ierr = PetscOptionsHead(PetscOptionsObject,"KSP FETIDP options");CHKERRQ(ierr); 1124 ierr = PetscOptionsBool("-ksp_fetidp_fullyredundant","Use fully redundant multipliers","none",fetidp->fully_redundant,&fetidp->fully_redundant,NULL);CHKERRQ(ierr); 1125 ierr = PetscOptionsBool("-ksp_fetidp_saddlepoint","Activates support for saddle-point problems",NULL,fetidp->saddlepoint,&fetidp->saddlepoint,NULL);CHKERRQ(ierr); 1126 ierr = PetscOptionsBool("-ksp_fetidp_check","Activates verbose debugging output FETI-DP operators",NULL,fetidp->check,&fetidp->check,NULL);CHKERRQ(ierr); 1127 ierr = PetscOptionsTail();CHKERRQ(ierr); 1128 ierr = PCSetFromOptions(fetidp->innerbddc);CHKERRQ(ierr); 1129 PetscFunctionReturn(0); 1130 } 1131 1132 /*MC 1133 KSPFETIDP - The FETI-DP method 1134 1135 This class implements the FETI-DP method [1]. 1136 The preconditioning matrix for the KSP must be of type MATIS. 1137 The FETI-DP linear system (automatically generated constructing an internal PCBDDC object) is solved using an inner KSP object. 1138 1139 Options Database Keys: 1140 + -ksp_fetidp_fullyredundant <false> : use a fully redundant set of Lagrange multipliers 1141 - -ksp_fetidp_saddlepoint <false> : activates support for saddle point problems, see [2] 1142 1143 Level: Advanced 1144 1145 Notes: Options for the inner KSP and for the customization of the PCBDDC object can be specified at command line by using the prefixes -fetidp_ and -fetidp_bddc_. E.g., 1146 .vb 1147 -fetidp_ksp_type gmres -fetidp_bddc_pc_bddc_symmetric false 1148 .ve 1149 will use GMRES for the solution of the linear system on the Lagrange multipliers, generated using a non-symmetric PCBDDC. 1150 For saddle point problems with continuous pressures, the operators for the interface pressure solver can be specified with KSPFETIDPSetPressureOperators(). 1151 1152 References: 1153 .vb 1154 . [1] - C. Farhat, M. Lesoinne, P. LeTallec, K. Pierson, and D. Rixen, FETI-DP: a dual-primal unified FETI method. I. A faster alternative to the two-level FETI method, Internat. J. Numer. Methods Engrg., 50 (2001), pp. 1523--1544 1155 . [2] - X. Tu, J. Li, A FETI-DP type domain decomposition algorithm for three-dimensional incompressible Stokes equations, SIAM J. Numer. Anal., 53 (2015), pp. 720-742 1156 .ve 1157 1158 .seealso: MATIS, PCBDDC, KSPFETIDPSetInnerBDDC, KSPFETIDPGetInnerBDDC, KSPFETIDPGetInnerKSP 1159 M*/ 1160 #undef __FUNCT__ 1161 #define __FUNCT__ "KSPCreate_FETIDP" 1162 PETSC_EXTERN PetscErrorCode KSPCreate_FETIDP(KSP ksp) 1163 { 1164 PetscErrorCode ierr; 1165 KSP_FETIDP *fetidp; 1166 KSP_FETIDPMon *monctx; 1167 PC_BDDC *pcbddc; 1168 PC pc; 1169 1170 PetscFunctionBegin; 1171 ierr = PetscNewLog(ksp,&fetidp);CHKERRQ(ierr); 1172 fetidp->matstate = -1; 1173 fetidp->matnnzstate = -1; 1174 fetidp->statechanged = PETSC_TRUE; 1175 1176 ksp->data = (void*)fetidp; 1177 ksp->ops->setup = KSPSetUp_FETIDP; 1178 ksp->ops->solve = KSPSolve_FETIDP; 1179 ksp->ops->destroy = KSPDestroy_FETIDP; 1180 ksp->ops->computeeigenvalues = KSPComputeEigenvalues_FETIDP; 1181 ksp->ops->computeextremesingularvalues = KSPComputeExtremeSingularValues_FETIDP; 1182 ksp->ops->view = KSPView_FETIDP; 1183 ksp->ops->setfromoptions = KSPSetFromOptions_FETIDP; 1184 ksp->ops->buildsolution = KSPBuildSolution_FETIDP; 1185 ksp->ops->buildresidual = KSPBuildResidualDefault; 1186 /* create the inner KSP for the Lagrange multipliers */ 1187 ierr = KSPCreate(PetscObjectComm((PetscObject)ksp),&fetidp->innerksp);CHKERRQ(ierr); 1188 ierr = KSPGetPC(fetidp->innerksp,&pc);CHKERRQ(ierr); 1189 ierr = PCSetType(pc,PCNONE);CHKERRQ(ierr); 1190 ierr = PetscLogObjectParent((PetscObject)ksp,(PetscObject)fetidp->innerksp);CHKERRQ(ierr); 1191 /* monitor */ 1192 ierr = PetscNew(&monctx);CHKERRQ(ierr); 1193 monctx->parentksp = ksp; 1194 fetidp->monctx = monctx; 1195 ierr = KSPMonitorSet(fetidp->innerksp,KSPMonitor_FETIDP,fetidp->monctx,NULL);CHKERRQ(ierr); 1196 /* create the inner BDDC */ 1197 ierr = PCCreate(PetscObjectComm((PetscObject)ksp),&fetidp->innerbddc);CHKERRQ(ierr); 1198 ierr = PCSetType(fetidp->innerbddc,PCBDDC);CHKERRQ(ierr); 1199 /* make sure we always obtain a consistent FETI-DP matrix 1200 for symmetric problems, the user can always customize it through the command line */ 1201 pcbddc = (PC_BDDC*)fetidp->innerbddc->data; 1202 pcbddc->symmetric_primal = PETSC_FALSE; 1203 ierr = PetscLogObjectParent((PetscObject)ksp,(PetscObject)fetidp->innerbddc);CHKERRQ(ierr); 1204 /* composed functions */ 1205 ierr = PetscObjectComposeFunction((PetscObject)ksp,"KSPFETIDPSetInnerBDDC_C",KSPFETIDPSetInnerBDDC_FETIDP);CHKERRQ(ierr); 1206 ierr = PetscObjectComposeFunction((PetscObject)ksp,"KSPFETIDPGetInnerBDDC_C",KSPFETIDPGetInnerBDDC_FETIDP);CHKERRQ(ierr); 1207 ierr = PetscObjectComposeFunction((PetscObject)ksp,"KSPFETIDPGetInnerKSP_C",KSPFETIDPGetInnerKSP_FETIDP);CHKERRQ(ierr); 1208 ierr = PetscObjectComposeFunction((PetscObject)ksp,"KSPFETIDPSetPressureOperators_C",KSPFETIDPSetPressureOperators_FETIDP);CHKERRQ(ierr); 1209 /* need to call KSPSetUp_FETIDP even with KSP_SETUP_NEWMATRIX */ 1210 ksp->setupnewmatrix = PETSC_TRUE; 1211 PetscFunctionReturn(0); 1212 } 1213