1 2 #include <../src/snes/impls/vi/rs/virsimpl.h> /*I "petscsnes.h" I*/ 3 #include <petsc/private/dmimpl.h> 4 #include <petsc/private/vecimpl.h> 5 6 /* 7 SNESVIGetInactiveSet - Gets the global indices for the inactive set variables (these correspond to the degrees of freedom the linear 8 system is solved on) 9 10 Input parameter: 11 . snes - the SNES context 12 13 Output parameter: 14 . inact - inactive set index set 15 16 */ 17 PetscErrorCode SNESVIGetInactiveSet(SNES snes,IS *inact) 18 { 19 SNES_VINEWTONRSLS *vi = (SNES_VINEWTONRSLS*)snes->data; 20 21 PetscFunctionBegin; 22 *inact = vi->IS_inact; 23 PetscFunctionReturn(0); 24 } 25 26 /* 27 Provides a wrapper to a DM to allow it to be used to generated the interpolation/restriction from the DM for the smaller matrices and vectors 28 defined by the reduced space method. 29 30 Simple calls the regular DM interpolation and restricts it to operation on the variables not associated with active constraints. 31 32 */ 33 typedef struct { 34 PetscInt n; /* size of vectors in the reduced DM space */ 35 IS inactive; 36 37 PetscErrorCode (*createinterpolation)(DM,DM,Mat*,Vec*); /* DM's original routines */ 38 PetscErrorCode (*coarsen)(DM, MPI_Comm, DM*); 39 PetscErrorCode (*createglobalvector)(DM,Vec*); 40 PetscErrorCode (*createinjection)(DM,DM,Mat*); 41 PetscErrorCode (*hascreateinjection)(DM,PetscBool*); 42 43 DM dm; /* when destroying this object we need to reset the above function into the base DM */ 44 } DM_SNESVI; 45 46 /* 47 DMCreateGlobalVector_SNESVI - Creates global vector of the size of the reduced space 48 49 */ 50 PetscErrorCode DMCreateGlobalVector_SNESVI(DM dm,Vec *vec) 51 { 52 PetscContainer isnes; 53 DM_SNESVI *dmsnesvi; 54 55 PetscFunctionBegin; 56 PetscCall(PetscObjectQuery((PetscObject)dm,"VI",(PetscObject*)&isnes)); 57 PetscCheck(isnes,PetscObjectComm((PetscObject)dm),PETSC_ERR_PLIB,"Composed SNES is missing"); 58 PetscCall(PetscContainerGetPointer(isnes,(void**)&dmsnesvi)); 59 PetscCall(VecCreateMPI(PetscObjectComm((PetscObject)dm),dmsnesvi->n,PETSC_DETERMINE,vec)); 60 PetscCall(VecSetDM(*vec, dm)); 61 PetscFunctionReturn(0); 62 } 63 64 /* 65 DMCreateInterpolation_SNESVI - Modifieds the interpolation obtained from the DM by removing all rows and columns associated with active constraints. 66 67 */ 68 PetscErrorCode DMCreateInterpolation_SNESVI(DM dm1,DM dm2,Mat *mat,Vec *vec) 69 { 70 PetscContainer isnes; 71 DM_SNESVI *dmsnesvi1,*dmsnesvi2; 72 Mat interp; 73 74 PetscFunctionBegin; 75 PetscCall(PetscObjectQuery((PetscObject)dm1,"VI",(PetscObject*)&isnes)); 76 PetscCheck(isnes,PetscObjectComm((PetscObject)dm1),PETSC_ERR_PLIB,"Composed VI data structure is missing"); 77 PetscCall(PetscContainerGetPointer(isnes,(void**)&dmsnesvi1)); 78 PetscCall(PetscObjectQuery((PetscObject)dm2,"VI",(PetscObject*)&isnes)); 79 PetscCheck(isnes,PetscObjectComm((PetscObject)dm2),PETSC_ERR_PLIB,"Composed VI data structure is missing"); 80 PetscCall(PetscContainerGetPointer(isnes,(void**)&dmsnesvi2)); 81 82 PetscCall((*dmsnesvi1->createinterpolation)(dm1,dm2,&interp,NULL)); 83 PetscCall(MatCreateSubMatrix(interp,dmsnesvi2->inactive,dmsnesvi1->inactive,MAT_INITIAL_MATRIX,mat)); 84 PetscCall(MatDestroy(&interp)); 85 *vec = NULL; 86 PetscFunctionReturn(0); 87 } 88 89 static PetscErrorCode DMSetVI(DM,IS); 90 static PetscErrorCode DMDestroyVI(DM); 91 92 /* 93 DMCoarsen_SNESVI - Computes the regular coarsened DM then computes additional information about its inactive set 94 95 */ 96 PetscErrorCode DMCoarsen_SNESVI(DM dm1,MPI_Comm comm,DM *dm2) 97 { 98 PetscContainer isnes; 99 DM_SNESVI *dmsnesvi1; 100 Vec finemarked,coarsemarked; 101 IS inactive; 102 Mat inject; 103 const PetscInt *index; 104 PetscInt n,k,cnt = 0,rstart,*coarseindex; 105 PetscScalar *marked; 106 107 PetscFunctionBegin; 108 PetscCall(PetscObjectQuery((PetscObject)dm1,"VI",(PetscObject*)&isnes)); 109 PetscCheck(isnes,PetscObjectComm((PetscObject)dm1),PETSC_ERR_PLIB,"Composed VI data structure is missing"); 110 PetscCall(PetscContainerGetPointer(isnes,(void**)&dmsnesvi1)); 111 112 /* get the original coarsen */ 113 PetscCall((*dmsnesvi1->coarsen)(dm1,comm,dm2)); 114 115 /* not sure why this extra reference is needed, but without the dm2 disappears too early */ 116 /* Updating the KSPCreateVecs() to avoid using DMGetGlobalVector() when matrix is available removes the need for this reference? */ 117 /* PetscCall(PetscObjectReference((PetscObject)*dm2));*/ 118 119 /* need to set back global vectors in order to use the original injection */ 120 PetscCall(DMClearGlobalVectors(dm1)); 121 122 dm1->ops->createglobalvector = dmsnesvi1->createglobalvector; 123 124 PetscCall(DMCreateGlobalVector(dm1,&finemarked)); 125 PetscCall(DMCreateGlobalVector(*dm2,&coarsemarked)); 126 127 /* 128 fill finemarked with locations of inactive points 129 */ 130 PetscCall(ISGetIndices(dmsnesvi1->inactive,&index)); 131 PetscCall(ISGetLocalSize(dmsnesvi1->inactive,&n)); 132 PetscCall(VecSet(finemarked,0.0)); 133 for (k=0; k<n; k++) { 134 PetscCall(VecSetValue(finemarked,index[k],1.0,INSERT_VALUES)); 135 } 136 PetscCall(VecAssemblyBegin(finemarked)); 137 PetscCall(VecAssemblyEnd(finemarked)); 138 139 PetscCall(DMCreateInjection(*dm2,dm1,&inject)); 140 PetscCall(MatRestrict(inject,finemarked,coarsemarked)); 141 PetscCall(MatDestroy(&inject)); 142 143 /* 144 create index set list of coarse inactive points from coarsemarked 145 */ 146 PetscCall(VecGetLocalSize(coarsemarked,&n)); 147 PetscCall(VecGetOwnershipRange(coarsemarked,&rstart,NULL)); 148 PetscCall(VecGetArray(coarsemarked,&marked)); 149 for (k=0; k<n; k++) { 150 if (marked[k] != 0.0) cnt++; 151 } 152 PetscCall(PetscMalloc1(cnt,&coarseindex)); 153 cnt = 0; 154 for (k=0; k<n; k++) { 155 if (marked[k] != 0.0) coarseindex[cnt++] = k + rstart; 156 } 157 PetscCall(VecRestoreArray(coarsemarked,&marked)); 158 PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)coarsemarked),cnt,coarseindex,PETSC_OWN_POINTER,&inactive)); 159 160 PetscCall(DMClearGlobalVectors(dm1)); 161 162 dm1->ops->createglobalvector = DMCreateGlobalVector_SNESVI; 163 164 PetscCall(DMSetVI(*dm2,inactive)); 165 166 PetscCall(VecDestroy(&finemarked)); 167 PetscCall(VecDestroy(&coarsemarked)); 168 PetscCall(ISDestroy(&inactive)); 169 PetscFunctionReturn(0); 170 } 171 172 PetscErrorCode DMDestroy_SNESVI(DM_SNESVI *dmsnesvi) 173 { 174 PetscFunctionBegin; 175 /* reset the base methods in the DM object that were changed when the DM_SNESVI was reset */ 176 dmsnesvi->dm->ops->createinterpolation = dmsnesvi->createinterpolation; 177 dmsnesvi->dm->ops->coarsen = dmsnesvi->coarsen; 178 dmsnesvi->dm->ops->createglobalvector = dmsnesvi->createglobalvector; 179 dmsnesvi->dm->ops->createinjection = dmsnesvi->createinjection; 180 dmsnesvi->dm->ops->hascreateinjection = dmsnesvi->hascreateinjection; 181 /* need to clear out this vectors because some of them may not have a reference to the DM 182 but they are counted as having references to the DM in DMDestroy() */ 183 PetscCall(DMClearGlobalVectors(dmsnesvi->dm)); 184 185 PetscCall(ISDestroy(&dmsnesvi->inactive)); 186 PetscCall(PetscFree(dmsnesvi)); 187 PetscFunctionReturn(0); 188 } 189 190 /* 191 DMSetVI - Marks a DM as associated with a VI problem. This causes the interpolation/restriction operators to 192 be restricted to only those variables NOT associated with active constraints. 193 194 */ 195 static PetscErrorCode DMSetVI(DM dm,IS inactive) 196 { 197 PetscContainer isnes; 198 DM_SNESVI *dmsnesvi; 199 200 PetscFunctionBegin; 201 if (!dm) PetscFunctionReturn(0); 202 203 PetscCall(PetscObjectReference((PetscObject)inactive)); 204 205 PetscCall(PetscObjectQuery((PetscObject)dm,"VI",(PetscObject*)&isnes)); 206 if (!isnes) { 207 PetscCall(PetscContainerCreate(PetscObjectComm((PetscObject)dm),&isnes)); 208 PetscCall(PetscContainerSetUserDestroy(isnes,(PetscErrorCode (*)(void*))DMDestroy_SNESVI)); 209 PetscCall(PetscNew(&dmsnesvi)); 210 PetscCall(PetscContainerSetPointer(isnes,(void*)dmsnesvi)); 211 PetscCall(PetscObjectCompose((PetscObject)dm,"VI",(PetscObject)isnes)); 212 PetscCall(PetscContainerDestroy(&isnes)); 213 214 dmsnesvi->createinterpolation = dm->ops->createinterpolation; 215 dm->ops->createinterpolation = DMCreateInterpolation_SNESVI; 216 dmsnesvi->coarsen = dm->ops->coarsen; 217 dm->ops->coarsen = DMCoarsen_SNESVI; 218 dmsnesvi->createglobalvector = dm->ops->createglobalvector; 219 dm->ops->createglobalvector = DMCreateGlobalVector_SNESVI; 220 dmsnesvi->createinjection = dm->ops->createinjection; 221 dm->ops->createinjection = NULL; 222 dmsnesvi->hascreateinjection = dm->ops->hascreateinjection; 223 dm->ops->hascreateinjection = NULL; 224 } else { 225 PetscCall(PetscContainerGetPointer(isnes,(void**)&dmsnesvi)); 226 PetscCall(ISDestroy(&dmsnesvi->inactive)); 227 } 228 PetscCall(DMClearGlobalVectors(dm)); 229 PetscCall(ISGetLocalSize(inactive,&dmsnesvi->n)); 230 231 dmsnesvi->inactive = inactive; 232 dmsnesvi->dm = dm; 233 PetscFunctionReturn(0); 234 } 235 236 /* 237 DMDestroyVI - Frees the DM_SNESVI object contained in the DM 238 - also resets the function pointers in the DM for createinterpolation() etc to use the original DM 239 */ 240 static PetscErrorCode DMDestroyVI(DM dm) 241 { 242 PetscFunctionBegin; 243 if (!dm) PetscFunctionReturn(0); 244 PetscCall(PetscObjectCompose((PetscObject)dm,"VI",(PetscObject)NULL)); 245 PetscFunctionReturn(0); 246 } 247 248 /* --------------------------------------------------------------------------------------------------------*/ 249 250 PetscErrorCode SNESCreateIndexSets_VINEWTONRSLS(SNES snes,Vec X,Vec F,IS *ISact,IS *ISinact) 251 { 252 PetscFunctionBegin; 253 PetscCall(SNESVIGetActiveSetIS(snes,X,F,ISact)); 254 PetscCall(ISComplement(*ISact,X->map->rstart,X->map->rend,ISinact)); 255 PetscFunctionReturn(0); 256 } 257 258 /* Create active and inactive set vectors. The local size of this vector is set and petsc computes the global size */ 259 PetscErrorCode SNESCreateSubVectors_VINEWTONRSLS(SNES snes,PetscInt n,Vec *newv) 260 { 261 Vec v; 262 263 PetscFunctionBegin; 264 PetscCall(VecCreate(PetscObjectComm((PetscObject)snes),&v)); 265 PetscCall(VecSetSizes(v,n,PETSC_DECIDE)); 266 PetscCall(VecSetType(v,VECSTANDARD)); 267 *newv = v; 268 PetscFunctionReturn(0); 269 } 270 271 /* Resets the snes PC and KSP when the active set sizes change */ 272 PetscErrorCode SNESVIResetPCandKSP(SNES snes,Mat Amat,Mat Pmat) 273 { 274 KSP snesksp; 275 276 PetscFunctionBegin; 277 PetscCall(SNESGetKSP(snes,&snesksp)); 278 PetscCall(KSPReset(snesksp)); 279 PetscCall(KSPResetFromOptions(snesksp)); 280 281 /* 282 KSP kspnew; 283 PC pcnew; 284 MatSolverType stype; 285 286 PetscCall(KSPCreate(PetscObjectComm((PetscObject)snes),&kspnew)); 287 kspnew->pc_side = snesksp->pc_side; 288 kspnew->rtol = snesksp->rtol; 289 kspnew->abstol = snesksp->abstol; 290 kspnew->max_it = snesksp->max_it; 291 PetscCall(KSPSetType(kspnew,((PetscObject)snesksp)->type_name)); 292 PetscCall(KSPGetPC(kspnew,&pcnew)); 293 PetscCall(PCSetType(kspnew->pc,((PetscObject)snesksp->pc)->type_name)); 294 PetscCall(PCSetOperators(kspnew->pc,Amat,Pmat)); 295 PetscCall(PCFactorGetMatSolverType(snesksp->pc,&stype)); 296 PetscCall(PCFactorSetMatSolverType(kspnew->pc,stype)); 297 PetscCall(KSPDestroy(&snesksp)); 298 snes->ksp = kspnew; 299 PetscCall(PetscLogObjectParent((PetscObject)snes,(PetscObject)kspnew)); 300 PetscCall(KSPSetFromOptions(kspnew));*/ 301 PetscFunctionReturn(0); 302 } 303 304 /* Variational Inequality solver using reduce space method. No semismooth algorithm is 305 implemented in this algorithm. It basically identifies the active constraints and does 306 a linear solve on the other variables (those not associated with the active constraints). */ 307 PetscErrorCode SNESSolve_VINEWTONRSLS(SNES snes) 308 { 309 SNES_VINEWTONRSLS *vi = (SNES_VINEWTONRSLS*)snes->data; 310 PetscInt maxits,i,lits; 311 SNESLineSearchReason lssucceed; 312 PetscReal fnorm,gnorm,xnorm=0,ynorm; 313 Vec Y,X,F; 314 KSPConvergedReason kspreason; 315 KSP ksp; 316 PC pc; 317 318 PetscFunctionBegin; 319 /* Multigrid must use Galerkin for coarse grids with active set/reduced space methods; cannot rediscretize on coarser grids*/ 320 PetscCall(SNESGetKSP(snes,&ksp)); 321 PetscCall(KSPGetPC(ksp,&pc)); 322 PetscCall(PCMGSetGalerkin(pc,PC_MG_GALERKIN_BOTH)); 323 324 snes->numFailures = 0; 325 snes->numLinearSolveFailures = 0; 326 snes->reason = SNES_CONVERGED_ITERATING; 327 328 maxits = snes->max_its; /* maximum number of iterations */ 329 X = snes->vec_sol; /* solution vector */ 330 F = snes->vec_func; /* residual vector */ 331 Y = snes->work[0]; /* work vectors */ 332 333 PetscCall(SNESLineSearchSetVIFunctions(snes->linesearch, SNESVIProjectOntoBounds, SNESVIComputeInactiveSetFnorm)); 334 PetscCall(SNESLineSearchSetVecs(snes->linesearch, X, NULL, NULL, NULL, NULL)); 335 PetscCall(SNESLineSearchSetUp(snes->linesearch)); 336 337 PetscCall(PetscObjectSAWsTakeAccess((PetscObject)snes)); 338 snes->iter = 0; 339 snes->norm = 0.0; 340 PetscCall(PetscObjectSAWsGrantAccess((PetscObject)snes)); 341 342 PetscCall(SNESVIProjectOntoBounds(snes,X)); 343 PetscCall(SNESComputeFunction(snes,X,F)); 344 PetscCall(SNESVIComputeInactiveSetFnorm(snes,F,X,&fnorm)); 345 PetscCall(VecNorm(X,NORM_2,&xnorm)); /* xnorm <- ||x|| */ 346 SNESCheckFunctionNorm(snes,fnorm); 347 PetscCall(PetscObjectSAWsTakeAccess((PetscObject)snes)); 348 snes->norm = fnorm; 349 PetscCall(PetscObjectSAWsGrantAccess((PetscObject)snes)); 350 PetscCall(SNESLogConvergenceHistory(snes,fnorm,0)); 351 PetscCall(SNESMonitor(snes,0,fnorm)); 352 353 /* test convergence */ 354 PetscCall((*snes->ops->converged)(snes,0,0.0,0.0,fnorm,&snes->reason,snes->cnvP)); 355 if (snes->reason) PetscFunctionReturn(0); 356 357 for (i=0; i<maxits; i++) { 358 359 IS IS_act; /* _act -> active set _inact -> inactive set */ 360 IS IS_redact; /* redundant active set */ 361 VecScatter scat_act,scat_inact; 362 PetscInt nis_act,nis_inact; 363 Vec Y_act,Y_inact,F_inact; 364 Mat jac_inact_inact,prejac_inact_inact; 365 PetscBool isequal; 366 367 /* Call general purpose update function */ 368 if (snes->ops->update) { 369 PetscCall((*snes->ops->update)(snes, snes->iter)); 370 } 371 PetscCall(SNESComputeJacobian(snes,X,snes->jacobian,snes->jacobian_pre)); 372 SNESCheckJacobianDomainerror(snes); 373 374 /* Create active and inactive index sets */ 375 376 /*original 377 PetscCall(SNESVICreateIndexSets_RS(snes,X,F,&IS_act,&vi->IS_inact)); 378 */ 379 PetscCall(SNESVIGetActiveSetIS(snes,X,F,&IS_act)); 380 381 if (vi->checkredundancy) { 382 PetscCall((*vi->checkredundancy)(snes,IS_act,&IS_redact,vi->ctxP)); 383 if (IS_redact) { 384 PetscCall(ISSort(IS_redact)); 385 PetscCall(ISComplement(IS_redact,X->map->rstart,X->map->rend,&vi->IS_inact)); 386 PetscCall(ISDestroy(&IS_redact)); 387 } else { 388 PetscCall(ISComplement(IS_act,X->map->rstart,X->map->rend,&vi->IS_inact)); 389 } 390 } else { 391 PetscCall(ISComplement(IS_act,X->map->rstart,X->map->rend,&vi->IS_inact)); 392 } 393 394 /* Create inactive set submatrix */ 395 PetscCall(MatCreateSubMatrix(snes->jacobian,vi->IS_inact,vi->IS_inact,MAT_INITIAL_MATRIX,&jac_inact_inact)); 396 397 if (0) { /* Dead code (temporary developer hack) */ 398 IS keptrows; 399 PetscCall(MatFindNonzeroRows(jac_inact_inact,&keptrows)); 400 if (keptrows) { 401 PetscInt cnt,*nrows,k; 402 const PetscInt *krows,*inact; 403 PetscInt rstart; 404 405 PetscCall(MatGetOwnershipRange(jac_inact_inact,&rstart,NULL)); 406 PetscCall(MatDestroy(&jac_inact_inact)); 407 PetscCall(ISDestroy(&IS_act)); 408 409 PetscCall(ISGetLocalSize(keptrows,&cnt)); 410 PetscCall(ISGetIndices(keptrows,&krows)); 411 PetscCall(ISGetIndices(vi->IS_inact,&inact)); 412 PetscCall(PetscMalloc1(cnt,&nrows)); 413 for (k=0; k<cnt; k++) nrows[k] = inact[krows[k]-rstart]; 414 PetscCall(ISRestoreIndices(keptrows,&krows)); 415 PetscCall(ISRestoreIndices(vi->IS_inact,&inact)); 416 PetscCall(ISDestroy(&keptrows)); 417 PetscCall(ISDestroy(&vi->IS_inact)); 418 419 PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)snes),cnt,nrows,PETSC_OWN_POINTER,&vi->IS_inact)); 420 PetscCall(ISComplement(vi->IS_inact,F->map->rstart,F->map->rend,&IS_act)); 421 PetscCall(MatCreateSubMatrix(snes->jacobian,vi->IS_inact,vi->IS_inact,MAT_INITIAL_MATRIX,&jac_inact_inact)); 422 } 423 } 424 PetscCall(DMSetVI(snes->dm,vi->IS_inact)); 425 /* remove later */ 426 427 /* 428 PetscCall(VecView(vi->xu,PETSC_VIEWER_BINARY_(((PetscObject)(vi->xu))->comm))); 429 PetscCall(VecView(vi->xl,PETSC_VIEWER_BINARY_(((PetscObject)(vi->xl))->comm))); 430 PetscCall(VecView(X,PETSC_VIEWER_BINARY_(PetscObjectComm((PetscObject)X)))); 431 PetscCall(VecView(F,PETSC_VIEWER_BINARY_(PetscObjectComm((PetscObject)F)))); 432 PetscCall(ISView(vi->IS_inact,PETSC_VIEWER_BINARY_(PetscObjectComm((PetscObject)vi->IS_inact)))); 433 */ 434 435 /* Get sizes of active and inactive sets */ 436 PetscCall(ISGetLocalSize(IS_act,&nis_act)); 437 PetscCall(ISGetLocalSize(vi->IS_inact,&nis_inact)); 438 439 /* Create active and inactive set vectors */ 440 PetscCall(SNESCreateSubVectors_VINEWTONRSLS(snes,nis_inact,&F_inact)); 441 PetscCall(SNESCreateSubVectors_VINEWTONRSLS(snes,nis_act,&Y_act)); 442 PetscCall(SNESCreateSubVectors_VINEWTONRSLS(snes,nis_inact,&Y_inact)); 443 444 /* Create scatter contexts */ 445 PetscCall(VecScatterCreate(Y,IS_act,Y_act,NULL,&scat_act)); 446 PetscCall(VecScatterCreate(Y,vi->IS_inact,Y_inact,NULL,&scat_inact)); 447 448 /* Do a vec scatter to active and inactive set vectors */ 449 PetscCall(VecScatterBegin(scat_inact,F,F_inact,INSERT_VALUES,SCATTER_FORWARD)); 450 PetscCall(VecScatterEnd(scat_inact,F,F_inact,INSERT_VALUES,SCATTER_FORWARD)); 451 452 PetscCall(VecScatterBegin(scat_act,Y,Y_act,INSERT_VALUES,SCATTER_FORWARD)); 453 PetscCall(VecScatterEnd(scat_act,Y,Y_act,INSERT_VALUES,SCATTER_FORWARD)); 454 455 PetscCall(VecScatterBegin(scat_inact,Y,Y_inact,INSERT_VALUES,SCATTER_FORWARD)); 456 PetscCall(VecScatterEnd(scat_inact,Y,Y_inact,INSERT_VALUES,SCATTER_FORWARD)); 457 458 /* Active set direction = 0 */ 459 PetscCall(VecSet(Y_act,0)); 460 if (snes->jacobian != snes->jacobian_pre) { 461 PetscCall(MatCreateSubMatrix(snes->jacobian_pre,vi->IS_inact,vi->IS_inact,MAT_INITIAL_MATRIX,&prejac_inact_inact)); 462 } else prejac_inact_inact = jac_inact_inact; 463 464 PetscCall(ISEqual(vi->IS_inact_prev,vi->IS_inact,&isequal)); 465 if (!isequal) { 466 PetscCall(SNESVIResetPCandKSP(snes,jac_inact_inact,prejac_inact_inact)); 467 PetscCall(PCFieldSplitRestrictIS(pc,vi->IS_inact)); 468 } 469 470 /* PetscCall(ISView(vi->IS_inact,0)); */ 471 /* PetscCall(ISView(IS_act,0));*/ 472 /* ierr = MatView(snes->jacobian_pre,0); */ 473 474 PetscCall(KSPSetOperators(snes->ksp,jac_inact_inact,prejac_inact_inact)); 475 PetscCall(KSPSetUp(snes->ksp)); 476 { 477 PC pc; 478 PetscBool flg; 479 PetscCall(KSPGetPC(snes->ksp,&pc)); 480 PetscCall(PetscObjectTypeCompare((PetscObject)pc,PCFIELDSPLIT,&flg)); 481 if (flg) { 482 KSP *subksps; 483 PetscCall(PCFieldSplitGetSubKSP(pc,NULL,&subksps)); 484 PetscCall(KSPGetPC(subksps[0],&pc)); 485 PetscCall(PetscFree(subksps)); 486 PetscCall(PetscObjectTypeCompare((PetscObject)pc,PCBJACOBI,&flg)); 487 if (flg) { 488 PetscInt n,N = 101*101,j,cnts[3] = {0,0,0}; 489 const PetscInt *ii; 490 491 PetscCall(ISGetSize(vi->IS_inact,&n)); 492 PetscCall(ISGetIndices(vi->IS_inact,&ii)); 493 for (j=0; j<n; j++) { 494 if (ii[j] < N) cnts[0]++; 495 else if (ii[j] < 2*N) cnts[1]++; 496 else if (ii[j] < 3*N) cnts[2]++; 497 } 498 PetscCall(ISRestoreIndices(vi->IS_inact,&ii)); 499 500 PetscCall(PCBJacobiSetTotalBlocks(pc,3,cnts)); 501 } 502 } 503 } 504 505 PetscCall(KSPSolve(snes->ksp,F_inact,Y_inact)); 506 PetscCall(VecScatterBegin(scat_act,Y_act,Y,INSERT_VALUES,SCATTER_REVERSE)); 507 PetscCall(VecScatterEnd(scat_act,Y_act,Y,INSERT_VALUES,SCATTER_REVERSE)); 508 PetscCall(VecScatterBegin(scat_inact,Y_inact,Y,INSERT_VALUES,SCATTER_REVERSE)); 509 PetscCall(VecScatterEnd(scat_inact,Y_inact,Y,INSERT_VALUES,SCATTER_REVERSE)); 510 511 PetscCall(VecDestroy(&F_inact)); 512 PetscCall(VecDestroy(&Y_act)); 513 PetscCall(VecDestroy(&Y_inact)); 514 PetscCall(VecScatterDestroy(&scat_act)); 515 PetscCall(VecScatterDestroy(&scat_inact)); 516 PetscCall(ISDestroy(&IS_act)); 517 if (!isequal) { 518 PetscCall(ISDestroy(&vi->IS_inact_prev)); 519 PetscCall(ISDuplicate(vi->IS_inact,&vi->IS_inact_prev)); 520 } 521 PetscCall(ISDestroy(&vi->IS_inact)); 522 PetscCall(MatDestroy(&jac_inact_inact)); 523 if (snes->jacobian != snes->jacobian_pre) { 524 PetscCall(MatDestroy(&prejac_inact_inact)); 525 } 526 527 PetscCall(KSPGetConvergedReason(snes->ksp,&kspreason)); 528 if (kspreason < 0) { 529 if (++snes->numLinearSolveFailures >= snes->maxLinearSolveFailures) { 530 PetscCall(PetscInfo(snes,"iter=%" PetscInt_FMT ", number linear solve failures %" PetscInt_FMT " greater than current SNES allowed, stopping solve\n",snes->iter,snes->numLinearSolveFailures)); 531 snes->reason = SNES_DIVERGED_LINEAR_SOLVE; 532 break; 533 } 534 } 535 536 PetscCall(KSPGetIterationNumber(snes->ksp,&lits)); 537 snes->linear_its += lits; 538 PetscCall(PetscInfo(snes,"iter=%" PetscInt_FMT ", linear solve iterations=%" PetscInt_FMT "\n",snes->iter,lits)); 539 /* 540 if (snes->ops->precheck) { 541 PetscBool changed_y = PETSC_FALSE; 542 PetscCall((*snes->ops->precheck)(snes,X,Y,snes->precheck,&changed_y)); 543 } 544 545 if (PetscLogPrintInfo) { 546 PetscCall(SNESVICheckResidual_Private(snes,snes->jacobian,F,Y,G,W)); 547 } 548 */ 549 /* Compute a (scaled) negative update in the line search routine: 550 Y <- X - lambda*Y 551 and evaluate G = function(Y) (depends on the line search). 552 */ 553 PetscCall(VecCopy(Y,snes->vec_sol_update)); 554 ynorm = 1; gnorm = fnorm; 555 PetscCall(SNESLineSearchApply(snes->linesearch, X, F, &gnorm, Y)); 556 PetscCall(SNESLineSearchGetReason(snes->linesearch, &lssucceed)); 557 PetscCall(SNESLineSearchGetNorms(snes->linesearch, &xnorm, &gnorm, &ynorm)); 558 PetscCall(PetscInfo(snes,"fnorm=%18.16e, gnorm=%18.16e, ynorm=%18.16e, lssucceed=%d\n",(double)fnorm,(double)gnorm,(double)ynorm,(int)lssucceed)); 559 if (snes->reason == SNES_DIVERGED_FUNCTION_COUNT) break; 560 if (snes->domainerror) { 561 snes->reason = SNES_DIVERGED_FUNCTION_DOMAIN; 562 PetscCall(DMDestroyVI(snes->dm)); 563 PetscFunctionReturn(0); 564 } 565 if (lssucceed) { 566 if (++snes->numFailures >= snes->maxFailures) { 567 PetscBool ismin; 568 snes->reason = SNES_DIVERGED_LINE_SEARCH; 569 PetscCall(SNESVICheckLocalMin_Private(snes,snes->jacobian,F,X,gnorm,&ismin)); 570 if (ismin) snes->reason = SNES_DIVERGED_LOCAL_MIN; 571 break; 572 } 573 } 574 PetscCall(DMDestroyVI(snes->dm)); 575 /* Update function and solution vectors */ 576 fnorm = gnorm; 577 /* Monitor convergence */ 578 PetscCall(PetscObjectSAWsTakeAccess((PetscObject)snes)); 579 snes->iter = i+1; 580 snes->norm = fnorm; 581 snes->xnorm = xnorm; 582 snes->ynorm = ynorm; 583 PetscCall(PetscObjectSAWsGrantAccess((PetscObject)snes)); 584 PetscCall(SNESLogConvergenceHistory(snes,snes->norm,lits)); 585 PetscCall(SNESMonitor(snes,snes->iter,snes->norm)); 586 /* Test for convergence, xnorm = || X || */ 587 if (snes->ops->converged != SNESConvergedSkip) PetscCall(VecNorm(X,NORM_2,&xnorm)); 588 PetscCall((*snes->ops->converged)(snes,snes->iter,xnorm,ynorm,fnorm,&snes->reason,snes->cnvP)); 589 if (snes->reason) break; 590 } 591 /* make sure that the VI information attached to the DM is removed if the for loop above was broken early due to some exceptional conditional */ 592 PetscCall(DMDestroyVI(snes->dm)); 593 if (i == maxits) { 594 PetscCall(PetscInfo(snes,"Maximum number of iterations has been reached: %" PetscInt_FMT "\n",maxits)); 595 if (!snes->reason) snes->reason = SNES_DIVERGED_MAX_IT; 596 } 597 PetscFunctionReturn(0); 598 } 599 600 PetscErrorCode SNESVISetRedundancyCheck(SNES snes,PetscErrorCode (*func)(SNES,IS,IS*,void*),void *ctx) 601 { 602 SNES_VINEWTONRSLS *vi = (SNES_VINEWTONRSLS*)snes->data; 603 604 PetscFunctionBegin; 605 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 606 vi->checkredundancy = func; 607 vi->ctxP = ctx; 608 PetscFunctionReturn(0); 609 } 610 611 #if defined(PETSC_HAVE_MATLAB_ENGINE) 612 #include <engine.h> 613 #include <mex.h> 614 typedef struct {char *funcname; mxArray *ctx;} SNESMatlabContext; 615 616 PetscErrorCode SNESVIRedundancyCheck_Matlab(SNES snes,IS is_act,IS *is_redact,void *ctx) 617 { 618 SNESMatlabContext *sctx = (SNESMatlabContext*)ctx; 619 int nlhs = 1, nrhs = 5; 620 mxArray *plhs[1], *prhs[5]; 621 long long int l1 = 0, l2 = 0, ls = 0; 622 PetscInt *indices=NULL; 623 624 PetscFunctionBegin; 625 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 626 PetscValidHeaderSpecific(is_act,IS_CLASSID,2); 627 PetscValidPointer(is_redact,3); 628 PetscCheckSameComm(snes,1,is_act,2); 629 630 /* Create IS for reduced active set of size 0, its size and indices will 631 bet set by the Matlab function */ 632 PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)snes),0,indices,PETSC_OWN_POINTER,is_redact)); 633 /* call Matlab function in ctx */ 634 PetscCall(PetscArraycpy(&ls,&snes,1)); 635 PetscCall(PetscArraycpy(&l1,&is_act,1)); 636 PetscCall(PetscArraycpy(&l2,is_redact,1)); 637 prhs[0] = mxCreateDoubleScalar((double)ls); 638 prhs[1] = mxCreateDoubleScalar((double)l1); 639 prhs[2] = mxCreateDoubleScalar((double)l2); 640 prhs[3] = mxCreateString(sctx->funcname); 641 prhs[4] = sctx->ctx; 642 PetscCall(mexCallMATLAB(nlhs,plhs,nrhs,prhs,"PetscSNESVIRedundancyCheckInternal")); 643 PetscCall(mxGetScalar(plhs[0])); 644 mxDestroyArray(prhs[0]); 645 mxDestroyArray(prhs[1]); 646 mxDestroyArray(prhs[2]); 647 mxDestroyArray(prhs[3]); 648 mxDestroyArray(plhs[0]); 649 PetscFunctionReturn(0); 650 } 651 652 PetscErrorCode SNESVISetRedundancyCheckMatlab(SNES snes,const char *func,mxArray *ctx) 653 { 654 SNESMatlabContext *sctx; 655 656 PetscFunctionBegin; 657 /* currently sctx is memory bleed */ 658 PetscCall(PetscNew(&sctx)); 659 PetscCall(PetscStrallocpy(func,&sctx->funcname)); 660 sctx->ctx = mxDuplicateArray(ctx); 661 PetscCall(SNESVISetRedundancyCheck(snes,SNESVIRedundancyCheck_Matlab,sctx)); 662 PetscFunctionReturn(0); 663 } 664 665 #endif 666 667 /* -------------------------------------------------------------------------- */ 668 /* 669 SNESSetUp_VINEWTONRSLS - Sets up the internal data structures for the later use 670 of the SNESVI nonlinear solver. 671 672 Input Parameter: 673 . snes - the SNES context 674 675 Application Interface Routine: SNESSetUp() 676 677 Notes: 678 For basic use of the SNES solvers, the user need not explicitly call 679 SNESSetUp(), since these actions will automatically occur during 680 the call to SNESSolve(). 681 */ 682 PetscErrorCode SNESSetUp_VINEWTONRSLS(SNES snes) 683 { 684 SNES_VINEWTONRSLS *vi = (SNES_VINEWTONRSLS*) snes->data; 685 PetscInt *indices; 686 PetscInt i,n,rstart,rend; 687 SNESLineSearch linesearch; 688 689 PetscFunctionBegin; 690 PetscCall(SNESSetUp_VI(snes)); 691 692 /* Set up previous active index set for the first snes solve 693 vi->IS_inact_prev = 0,1,2,....N */ 694 695 PetscCall(VecGetOwnershipRange(snes->vec_sol,&rstart,&rend)); 696 PetscCall(VecGetLocalSize(snes->vec_sol,&n)); 697 PetscCall(PetscMalloc1(n,&indices)); 698 for (i=0; i < n; i++) indices[i] = rstart + i; 699 PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)snes),n,indices,PETSC_OWN_POINTER,&vi->IS_inact_prev)); 700 701 /* set the line search functions */ 702 if (!snes->linesearch) { 703 PetscCall(SNESGetLineSearch(snes, &linesearch)); 704 PetscCall(SNESLineSearchSetType(linesearch, SNESLINESEARCHBT)); 705 } 706 PetscFunctionReturn(0); 707 } 708 /* -------------------------------------------------------------------------- */ 709 PetscErrorCode SNESReset_VINEWTONRSLS(SNES snes) 710 { 711 SNES_VINEWTONRSLS *vi = (SNES_VINEWTONRSLS*) snes->data; 712 713 PetscFunctionBegin; 714 PetscCall(SNESReset_VI(snes)); 715 PetscCall(ISDestroy(&vi->IS_inact_prev)); 716 PetscFunctionReturn(0); 717 } 718 719 /* -------------------------------------------------------------------------- */ 720 /*MC 721 SNESVINEWTONRSLS - Reduced space active set solvers for variational inequalities based on Newton's method 722 723 Options Database: 724 + -snes_type <vinewtonssls,vinewtonrsls> - a semi-smooth solver, a reduced space active set method 725 - -snes_vi_monitor - prints the number of active constraints at each iteration. 726 727 Level: beginner 728 729 References: 730 . * - T. S. Munson, and S. Benson. Flexible Complementarity Solvers for Large Scale 731 Applications, Optimization Methods and Software, 21 (2006). 732 733 .seealso: `SNESVISetVariableBounds()`, `SNESVISetComputeVariableBounds()`, `SNESCreate()`, `SNES`, `SNESSetType()`, `SNESVINEWTONSSLS`, `SNESNEWTONTR`, `SNESLineSearchSetType()`, `SNESLineSearchSetPostCheck()`, `SNESLineSearchSetPreCheck()` 734 735 M*/ 736 PETSC_EXTERN PetscErrorCode SNESCreate_VINEWTONRSLS(SNES snes) 737 { 738 SNES_VINEWTONRSLS *vi; 739 SNESLineSearch linesearch; 740 741 PetscFunctionBegin; 742 snes->ops->reset = SNESReset_VINEWTONRSLS; 743 snes->ops->setup = SNESSetUp_VINEWTONRSLS; 744 snes->ops->solve = SNESSolve_VINEWTONRSLS; 745 snes->ops->destroy = SNESDestroy_VI; 746 snes->ops->setfromoptions = SNESSetFromOptions_VI; 747 snes->ops->view = NULL; 748 snes->ops->converged = SNESConvergedDefault_VI; 749 750 snes->usesksp = PETSC_TRUE; 751 snes->usesnpc = PETSC_FALSE; 752 753 PetscCall(SNESGetLineSearch(snes, &linesearch)); 754 if (!((PetscObject)linesearch)->type_name) { 755 PetscCall(SNESLineSearchSetType(linesearch, SNESLINESEARCHBT)); 756 } 757 PetscCall(SNESLineSearchBTSetAlpha(linesearch, 0.0)); 758 759 snes->alwayscomputesfinalresidual = PETSC_TRUE; 760 761 PetscCall(PetscNewLog(snes,&vi)); 762 snes->data = (void*)vi; 763 vi->checkredundancy = NULL; 764 765 PetscCall(PetscObjectComposeFunction((PetscObject)snes,"SNESVISetVariableBounds_C",SNESVISetVariableBounds_VI)); 766 PetscCall(PetscObjectComposeFunction((PetscObject)snes,"SNESVISetComputeVariableBounds_C",SNESVISetComputeVariableBounds_VI)); 767 PetscFunctionReturn(0); 768 } 769