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