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) PetscCall((*snes->ops->update)(snes, snes->iter)); 369 PetscCall(SNESComputeJacobian(snes,X,snes->jacobian,snes->jacobian_pre)); 370 SNESCheckJacobianDomainerror(snes); 371 372 /* Create active and inactive index sets */ 373 374 /*original 375 PetscCall(SNESVICreateIndexSets_RS(snes,X,F,&IS_act,&vi->IS_inact)); 376 */ 377 PetscCall(SNESVIGetActiveSetIS(snes,X,F,&IS_act)); 378 379 if (vi->checkredundancy) { 380 PetscCall((*vi->checkredundancy)(snes,IS_act,&IS_redact,vi->ctxP)); 381 if (IS_redact) { 382 PetscCall(ISSort(IS_redact)); 383 PetscCall(ISComplement(IS_redact,X->map->rstart,X->map->rend,&vi->IS_inact)); 384 PetscCall(ISDestroy(&IS_redact)); 385 } else { 386 PetscCall(ISComplement(IS_act,X->map->rstart,X->map->rend,&vi->IS_inact)); 387 } 388 } else { 389 PetscCall(ISComplement(IS_act,X->map->rstart,X->map->rend,&vi->IS_inact)); 390 } 391 392 /* Create inactive set submatrix */ 393 PetscCall(MatCreateSubMatrix(snes->jacobian,vi->IS_inact,vi->IS_inact,MAT_INITIAL_MATRIX,&jac_inact_inact)); 394 395 if (0) { /* Dead code (temporary developer hack) */ 396 IS keptrows; 397 PetscCall(MatFindNonzeroRows(jac_inact_inact,&keptrows)); 398 if (keptrows) { 399 PetscInt cnt,*nrows,k; 400 const PetscInt *krows,*inact; 401 PetscInt rstart; 402 403 PetscCall(MatGetOwnershipRange(jac_inact_inact,&rstart,NULL)); 404 PetscCall(MatDestroy(&jac_inact_inact)); 405 PetscCall(ISDestroy(&IS_act)); 406 407 PetscCall(ISGetLocalSize(keptrows,&cnt)); 408 PetscCall(ISGetIndices(keptrows,&krows)); 409 PetscCall(ISGetIndices(vi->IS_inact,&inact)); 410 PetscCall(PetscMalloc1(cnt,&nrows)); 411 for (k=0; k<cnt; k++) nrows[k] = inact[krows[k]-rstart]; 412 PetscCall(ISRestoreIndices(keptrows,&krows)); 413 PetscCall(ISRestoreIndices(vi->IS_inact,&inact)); 414 PetscCall(ISDestroy(&keptrows)); 415 PetscCall(ISDestroy(&vi->IS_inact)); 416 417 PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)snes),cnt,nrows,PETSC_OWN_POINTER,&vi->IS_inact)); 418 PetscCall(ISComplement(vi->IS_inact,F->map->rstart,F->map->rend,&IS_act)); 419 PetscCall(MatCreateSubMatrix(snes->jacobian,vi->IS_inact,vi->IS_inact,MAT_INITIAL_MATRIX,&jac_inact_inact)); 420 } 421 } 422 PetscCall(DMSetVI(snes->dm,vi->IS_inact)); 423 /* remove later */ 424 425 /* 426 PetscCall(VecView(vi->xu,PETSC_VIEWER_BINARY_(((PetscObject)(vi->xu))->comm))); 427 PetscCall(VecView(vi->xl,PETSC_VIEWER_BINARY_(((PetscObject)(vi->xl))->comm))); 428 PetscCall(VecView(X,PETSC_VIEWER_BINARY_(PetscObjectComm((PetscObject)X)))); 429 PetscCall(VecView(F,PETSC_VIEWER_BINARY_(PetscObjectComm((PetscObject)F)))); 430 PetscCall(ISView(vi->IS_inact,PETSC_VIEWER_BINARY_(PetscObjectComm((PetscObject)vi->IS_inact)))); 431 */ 432 433 /* Get sizes of active and inactive sets */ 434 PetscCall(ISGetLocalSize(IS_act,&nis_act)); 435 PetscCall(ISGetLocalSize(vi->IS_inact,&nis_inact)); 436 437 /* Create active and inactive set vectors */ 438 PetscCall(SNESCreateSubVectors_VINEWTONRSLS(snes,nis_inact,&F_inact)); 439 PetscCall(SNESCreateSubVectors_VINEWTONRSLS(snes,nis_act,&Y_act)); 440 PetscCall(SNESCreateSubVectors_VINEWTONRSLS(snes,nis_inact,&Y_inact)); 441 442 /* Create scatter contexts */ 443 PetscCall(VecScatterCreate(Y,IS_act,Y_act,NULL,&scat_act)); 444 PetscCall(VecScatterCreate(Y,vi->IS_inact,Y_inact,NULL,&scat_inact)); 445 446 /* Do a vec scatter to active and inactive set vectors */ 447 PetscCall(VecScatterBegin(scat_inact,F,F_inact,INSERT_VALUES,SCATTER_FORWARD)); 448 PetscCall(VecScatterEnd(scat_inact,F,F_inact,INSERT_VALUES,SCATTER_FORWARD)); 449 450 PetscCall(VecScatterBegin(scat_act,Y,Y_act,INSERT_VALUES,SCATTER_FORWARD)); 451 PetscCall(VecScatterEnd(scat_act,Y,Y_act,INSERT_VALUES,SCATTER_FORWARD)); 452 453 PetscCall(VecScatterBegin(scat_inact,Y,Y_inact,INSERT_VALUES,SCATTER_FORWARD)); 454 PetscCall(VecScatterEnd(scat_inact,Y,Y_inact,INSERT_VALUES,SCATTER_FORWARD)); 455 456 /* Active set direction = 0 */ 457 PetscCall(VecSet(Y_act,0)); 458 if (snes->jacobian != snes->jacobian_pre) { 459 PetscCall(MatCreateSubMatrix(snes->jacobian_pre,vi->IS_inact,vi->IS_inact,MAT_INITIAL_MATRIX,&prejac_inact_inact)); 460 } else prejac_inact_inact = jac_inact_inact; 461 462 PetscCall(ISEqual(vi->IS_inact_prev,vi->IS_inact,&isequal)); 463 if (!isequal) { 464 PetscCall(SNESVIResetPCandKSP(snes,jac_inact_inact,prejac_inact_inact)); 465 PetscCall(PCFieldSplitRestrictIS(pc,vi->IS_inact)); 466 } 467 468 /* PetscCall(ISView(vi->IS_inact,0)); */ 469 /* PetscCall(ISView(IS_act,0));*/ 470 /* ierr = MatView(snes->jacobian_pre,0); */ 471 472 PetscCall(KSPSetOperators(snes->ksp,jac_inact_inact,prejac_inact_inact)); 473 PetscCall(KSPSetUp(snes->ksp)); 474 { 475 PC pc; 476 PetscBool flg; 477 PetscCall(KSPGetPC(snes->ksp,&pc)); 478 PetscCall(PetscObjectTypeCompare((PetscObject)pc,PCFIELDSPLIT,&flg)); 479 if (flg) { 480 KSP *subksps; 481 PetscCall(PCFieldSplitGetSubKSP(pc,NULL,&subksps)); 482 PetscCall(KSPGetPC(subksps[0],&pc)); 483 PetscCall(PetscFree(subksps)); 484 PetscCall(PetscObjectTypeCompare((PetscObject)pc,PCBJACOBI,&flg)); 485 if (flg) { 486 PetscInt n,N = 101*101,j,cnts[3] = {0,0,0}; 487 const PetscInt *ii; 488 489 PetscCall(ISGetSize(vi->IS_inact,&n)); 490 PetscCall(ISGetIndices(vi->IS_inact,&ii)); 491 for (j=0; j<n; j++) { 492 if (ii[j] < N) cnts[0]++; 493 else if (ii[j] < 2*N) cnts[1]++; 494 else if (ii[j] < 3*N) cnts[2]++; 495 } 496 PetscCall(ISRestoreIndices(vi->IS_inact,&ii)); 497 498 PetscCall(PCBJacobiSetTotalBlocks(pc,3,cnts)); 499 } 500 } 501 } 502 503 PetscCall(KSPSolve(snes->ksp,F_inact,Y_inact)); 504 PetscCall(VecScatterBegin(scat_act,Y_act,Y,INSERT_VALUES,SCATTER_REVERSE)); 505 PetscCall(VecScatterEnd(scat_act,Y_act,Y,INSERT_VALUES,SCATTER_REVERSE)); 506 PetscCall(VecScatterBegin(scat_inact,Y_inact,Y,INSERT_VALUES,SCATTER_REVERSE)); 507 PetscCall(VecScatterEnd(scat_inact,Y_inact,Y,INSERT_VALUES,SCATTER_REVERSE)); 508 509 PetscCall(VecDestroy(&F_inact)); 510 PetscCall(VecDestroy(&Y_act)); 511 PetscCall(VecDestroy(&Y_inact)); 512 PetscCall(VecScatterDestroy(&scat_act)); 513 PetscCall(VecScatterDestroy(&scat_inact)); 514 PetscCall(ISDestroy(&IS_act)); 515 if (!isequal) { 516 PetscCall(ISDestroy(&vi->IS_inact_prev)); 517 PetscCall(ISDuplicate(vi->IS_inact,&vi->IS_inact_prev)); 518 } 519 PetscCall(ISDestroy(&vi->IS_inact)); 520 PetscCall(MatDestroy(&jac_inact_inact)); 521 if (snes->jacobian != snes->jacobian_pre) { 522 PetscCall(MatDestroy(&prejac_inact_inact)); 523 } 524 525 PetscCall(KSPGetConvergedReason(snes->ksp,&kspreason)); 526 if (kspreason < 0) { 527 if (++snes->numLinearSolveFailures >= snes->maxLinearSolveFailures) { 528 PetscCall(PetscInfo(snes,"iter=%" PetscInt_FMT ", number linear solve failures %" PetscInt_FMT " greater than current SNES allowed, stopping solve\n",snes->iter,snes->numLinearSolveFailures)); 529 snes->reason = SNES_DIVERGED_LINEAR_SOLVE; 530 break; 531 } 532 } 533 534 PetscCall(KSPGetIterationNumber(snes->ksp,&lits)); 535 snes->linear_its += lits; 536 PetscCall(PetscInfo(snes,"iter=%" PetscInt_FMT ", linear solve iterations=%" PetscInt_FMT "\n",snes->iter,lits)); 537 /* 538 if (snes->ops->precheck) { 539 PetscBool changed_y = PETSC_FALSE; 540 PetscCall((*snes->ops->precheck)(snes,X,Y,snes->precheck,&changed_y)); 541 } 542 543 if (PetscLogPrintInfo) PetscCall(SNESVICheckResidual_Private(snes,snes->jacobian,F,Y,G,W)); 544 */ 545 /* Compute a (scaled) negative update in the line search routine: 546 Y <- X - lambda*Y 547 and evaluate G = function(Y) (depends on the line search). 548 */ 549 PetscCall(VecCopy(Y,snes->vec_sol_update)); 550 ynorm = 1; gnorm = fnorm; 551 PetscCall(SNESLineSearchApply(snes->linesearch, X, F, &gnorm, Y)); 552 PetscCall(SNESLineSearchGetReason(snes->linesearch, &lssucceed)); 553 PetscCall(SNESLineSearchGetNorms(snes->linesearch, &xnorm, &gnorm, &ynorm)); 554 PetscCall(PetscInfo(snes,"fnorm=%18.16e, gnorm=%18.16e, ynorm=%18.16e, lssucceed=%d\n",(double)fnorm,(double)gnorm,(double)ynorm,(int)lssucceed)); 555 if (snes->reason == SNES_DIVERGED_FUNCTION_COUNT) break; 556 if (snes->domainerror) { 557 snes->reason = SNES_DIVERGED_FUNCTION_DOMAIN; 558 PetscCall(DMDestroyVI(snes->dm)); 559 PetscFunctionReturn(0); 560 } 561 if (lssucceed) { 562 if (++snes->numFailures >= snes->maxFailures) { 563 PetscBool ismin; 564 snes->reason = SNES_DIVERGED_LINE_SEARCH; 565 PetscCall(SNESVICheckLocalMin_Private(snes,snes->jacobian,F,X,gnorm,&ismin)); 566 if (ismin) snes->reason = SNES_DIVERGED_LOCAL_MIN; 567 break; 568 } 569 } 570 PetscCall(DMDestroyVI(snes->dm)); 571 /* Update function and solution vectors */ 572 fnorm = gnorm; 573 /* Monitor convergence */ 574 PetscCall(PetscObjectSAWsTakeAccess((PetscObject)snes)); 575 snes->iter = i+1; 576 snes->norm = fnorm; 577 snes->xnorm = xnorm; 578 snes->ynorm = ynorm; 579 PetscCall(PetscObjectSAWsGrantAccess((PetscObject)snes)); 580 PetscCall(SNESLogConvergenceHistory(snes,snes->norm,lits)); 581 PetscCall(SNESMonitor(snes,snes->iter,snes->norm)); 582 /* Test for convergence, xnorm = || X || */ 583 if (snes->ops->converged != SNESConvergedSkip) PetscCall(VecNorm(X,NORM_2,&xnorm)); 584 PetscCall((*snes->ops->converged)(snes,snes->iter,xnorm,ynorm,fnorm,&snes->reason,snes->cnvP)); 585 if (snes->reason) break; 586 } 587 /* 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 */ 588 PetscCall(DMDestroyVI(snes->dm)); 589 if (i == maxits) { 590 PetscCall(PetscInfo(snes,"Maximum number of iterations has been reached: %" PetscInt_FMT "\n",maxits)); 591 if (!snes->reason) snes->reason = SNES_DIVERGED_MAX_IT; 592 } 593 PetscFunctionReturn(0); 594 } 595 596 PetscErrorCode SNESVISetRedundancyCheck(SNES snes,PetscErrorCode (*func)(SNES,IS,IS*,void*),void *ctx) 597 { 598 SNES_VINEWTONRSLS *vi = (SNES_VINEWTONRSLS*)snes->data; 599 600 PetscFunctionBegin; 601 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 602 vi->checkredundancy = func; 603 vi->ctxP = ctx; 604 PetscFunctionReturn(0); 605 } 606 607 #if defined(PETSC_HAVE_MATLAB_ENGINE) 608 #include <engine.h> 609 #include <mex.h> 610 typedef struct {char *funcname; mxArray *ctx;} SNESMatlabContext; 611 612 PetscErrorCode SNESVIRedundancyCheck_Matlab(SNES snes,IS is_act,IS *is_redact,void *ctx) 613 { 614 SNESMatlabContext *sctx = (SNESMatlabContext*)ctx; 615 int nlhs = 1, nrhs = 5; 616 mxArray *plhs[1], *prhs[5]; 617 long long int l1 = 0, l2 = 0, ls = 0; 618 PetscInt *indices=NULL; 619 620 PetscFunctionBegin; 621 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 622 PetscValidHeaderSpecific(is_act,IS_CLASSID,2); 623 PetscValidPointer(is_redact,3); 624 PetscCheckSameComm(snes,1,is_act,2); 625 626 /* Create IS for reduced active set of size 0, its size and indices will 627 bet set by the Matlab function */ 628 PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)snes),0,indices,PETSC_OWN_POINTER,is_redact)); 629 /* call Matlab function in ctx */ 630 PetscCall(PetscArraycpy(&ls,&snes,1)); 631 PetscCall(PetscArraycpy(&l1,&is_act,1)); 632 PetscCall(PetscArraycpy(&l2,is_redact,1)); 633 prhs[0] = mxCreateDoubleScalar((double)ls); 634 prhs[1] = mxCreateDoubleScalar((double)l1); 635 prhs[2] = mxCreateDoubleScalar((double)l2); 636 prhs[3] = mxCreateString(sctx->funcname); 637 prhs[4] = sctx->ctx; 638 PetscCall(mexCallMATLAB(nlhs,plhs,nrhs,prhs,"PetscSNESVIRedundancyCheckInternal")); 639 PetscCall(mxGetScalar(plhs[0])); 640 mxDestroyArray(prhs[0]); 641 mxDestroyArray(prhs[1]); 642 mxDestroyArray(prhs[2]); 643 mxDestroyArray(prhs[3]); 644 mxDestroyArray(plhs[0]); 645 PetscFunctionReturn(0); 646 } 647 648 PetscErrorCode SNESVISetRedundancyCheckMatlab(SNES snes,const char *func,mxArray *ctx) 649 { 650 SNESMatlabContext *sctx; 651 652 PetscFunctionBegin; 653 /* currently sctx is memory bleed */ 654 PetscCall(PetscNew(&sctx)); 655 PetscCall(PetscStrallocpy(func,&sctx->funcname)); 656 sctx->ctx = mxDuplicateArray(ctx); 657 PetscCall(SNESVISetRedundancyCheck(snes,SNESVIRedundancyCheck_Matlab,sctx)); 658 PetscFunctionReturn(0); 659 } 660 661 #endif 662 663 /* -------------------------------------------------------------------------- */ 664 /* 665 SNESSetUp_VINEWTONRSLS - Sets up the internal data structures for the later use 666 of the SNESVI nonlinear solver. 667 668 Input Parameter: 669 . snes - the SNES context 670 671 Application Interface Routine: SNESSetUp() 672 673 Notes: 674 For basic use of the SNES solvers, the user need not explicitly call 675 SNESSetUp(), since these actions will automatically occur during 676 the call to SNESSolve(). 677 */ 678 PetscErrorCode SNESSetUp_VINEWTONRSLS(SNES snes) 679 { 680 SNES_VINEWTONRSLS *vi = (SNES_VINEWTONRSLS*) snes->data; 681 PetscInt *indices; 682 PetscInt i,n,rstart,rend; 683 SNESLineSearch linesearch; 684 685 PetscFunctionBegin; 686 PetscCall(SNESSetUp_VI(snes)); 687 688 /* Set up previous active index set for the first snes solve 689 vi->IS_inact_prev = 0,1,2,....N */ 690 691 PetscCall(VecGetOwnershipRange(snes->vec_sol,&rstart,&rend)); 692 PetscCall(VecGetLocalSize(snes->vec_sol,&n)); 693 PetscCall(PetscMalloc1(n,&indices)); 694 for (i=0; i < n; i++) indices[i] = rstart + i; 695 PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)snes),n,indices,PETSC_OWN_POINTER,&vi->IS_inact_prev)); 696 697 /* set the line search functions */ 698 if (!snes->linesearch) { 699 PetscCall(SNESGetLineSearch(snes, &linesearch)); 700 PetscCall(SNESLineSearchSetType(linesearch, SNESLINESEARCHBT)); 701 } 702 PetscFunctionReturn(0); 703 } 704 /* -------------------------------------------------------------------------- */ 705 PetscErrorCode SNESReset_VINEWTONRSLS(SNES snes) 706 { 707 SNES_VINEWTONRSLS *vi = (SNES_VINEWTONRSLS*) snes->data; 708 709 PetscFunctionBegin; 710 PetscCall(SNESReset_VI(snes)); 711 PetscCall(ISDestroy(&vi->IS_inact_prev)); 712 PetscFunctionReturn(0); 713 } 714 715 /* -------------------------------------------------------------------------- */ 716 /*MC 717 SNESVINEWTONRSLS - Reduced space active set solvers for variational inequalities based on Newton's method 718 719 Options Database: 720 + -snes_type <vinewtonssls,vinewtonrsls> - a semi-smooth solver, a reduced space active set method 721 - -snes_vi_monitor - prints the number of active constraints at each iteration. 722 723 Level: beginner 724 725 References: 726 . * - T. S. Munson, and S. Benson. Flexible Complementarity Solvers for Large Scale 727 Applications, Optimization Methods and Software, 21 (2006). 728 729 .seealso: `SNESVISetVariableBounds()`, `SNESVISetComputeVariableBounds()`, `SNESCreate()`, `SNES`, `SNESSetType()`, `SNESVINEWTONSSLS`, `SNESNEWTONTR`, `SNESLineSearchSetType()`, `SNESLineSearchSetPostCheck()`, `SNESLineSearchSetPreCheck()` 730 731 M*/ 732 PETSC_EXTERN PetscErrorCode SNESCreate_VINEWTONRSLS(SNES snes) 733 { 734 SNES_VINEWTONRSLS *vi; 735 SNESLineSearch linesearch; 736 737 PetscFunctionBegin; 738 snes->ops->reset = SNESReset_VINEWTONRSLS; 739 snes->ops->setup = SNESSetUp_VINEWTONRSLS; 740 snes->ops->solve = SNESSolve_VINEWTONRSLS; 741 snes->ops->destroy = SNESDestroy_VI; 742 snes->ops->setfromoptions = SNESSetFromOptions_VI; 743 snes->ops->view = NULL; 744 snes->ops->converged = SNESConvergedDefault_VI; 745 746 snes->usesksp = PETSC_TRUE; 747 snes->usesnpc = PETSC_FALSE; 748 749 PetscCall(SNESGetLineSearch(snes, &linesearch)); 750 if (!((PetscObject)linesearch)->type_name) { 751 PetscCall(SNESLineSearchSetType(linesearch, SNESLINESEARCHBT)); 752 } 753 PetscCall(SNESLineSearchBTSetAlpha(linesearch, 0.0)); 754 755 snes->alwayscomputesfinalresidual = PETSC_TRUE; 756 757 PetscCall(PetscNewLog(snes,&vi)); 758 snes->data = (void*)vi; 759 vi->checkredundancy = NULL; 760 761 PetscCall(PetscObjectComposeFunction((PetscObject)snes,"SNESVISetVariableBounds_C",SNESVISetVariableBounds_VI)); 762 PetscCall(PetscObjectComposeFunction((PetscObject)snes,"SNESVISetComputeVariableBounds_C",SNESVISetComputeVariableBounds_VI)); 763 PetscFunctionReturn(0); 764 } 765