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